PhreeqcRM coupling with a transport model

[hamiderfani]

Dear All,
I have a transport code and want to couple PhreeqcRM with my code to take account of rock-fluid interactions. I have previously coupled Iphreeqc successfully and it works fine but the runtime is high so I decided to replace it with PhreeqcRM.
I have different rock types (clay, quartz, calcite etc. ) in a model and some blocks are kinetic and some others are equilibrium controlled. As there is no clear and thorough manual and examples for PhreeqcRM I am a little bit confused. I just want to update my concentrations with PhreeqcRM and do not need its transport capabilities and my reaction spots are too much that I cannot define them by hand. I want to insert concentration and reactive surface area and time step at each block and receive updated ion concentrations based on the rock type of the block in addition to amount of dissolved or precipitated mineral. Can any body kindly describe the algorithm and needed functions for me?
should I define different datafiles based on the number of rock types? (one datafile for kinetic clay reaction with the rate expression and another for equilibrium calcite and  ...)
how can I insert the ion concentrations and the reactive surface area of mineral at each data block?

Thank you in advance for your response.
Hamid

[dlparkhurst]

Note that PhreeqcRM uses IPhreeqc to do the calculations, so it is a matter of using IPhreeqc efficiently. PhreeqcRM provides methods that are commonly used for reactive-transport calculation and includes parallelization.

Documentation of all methods is distributed with PhreeqcRM and is also available at

https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqcrm/index.html

An advection example in C, Fortran, and C++ is distributed with the code and demonstrates the use of most of the methods.

All of the capabilities you describe are available, documented, and shown in the example, except updating surface area. To do that, you probably need to resort to SURFACE_MODIFY blocks. Some care will be needed to generate a set of SURFACE_MODIFY blocks for the set of cells in each of the parallel threads (OpenMP) or processes (MPI).

[hamiderfani]

Dear Dr. Parkhurst
Thank you so much for the reponse, for the time being I am studying the privided example with the distribution of PhreeqcRM to learn how to use it for my own case but there are 2 problems
1. When I see the concentrations at the end of "set initial condition section" where the system is supposed to be initialized by solution 1, H+ concentration is 2.075714746752055E-007 despite pH is set to 7 in solution 1. Despite the fact that the concentration unit is set as mol/L (option 2) when I change the representative volume for the cells the concentration of H+ (and other species) changes in concentration matrix.
2. when I decrese the representative volume of the cells to 0.1 all concentration obtained by RM_GetConcentrations(id, c) are multiplied by a factor of 10 which is strange to me as the concentration must not change respect to representative volume and still no reaction is happenening in the system.

[dlparkhurst]

Using RV=1, I get pH 6.997. So, I don't know quite where you are getting your number. (I used SetPrintChemistryOn(true,false,false) before the RunCells in the initialization section to send results to advect_cpp.chem.txt.)

The second is unfortunately a bug. InitialPhreeqc2Module does not multiply by the rv factor as it should. So, until we put out a new release, you should use this method only with rv=1.0. InitialPhreeqcCell2Module works correctly with nonunity rv.

[hamiderfani]

Dear Dr. Parkhurst Thank you so much for the response,
As I am modeling multicomponent diffusive transport I understood that it is better to retrieve and import concentrations using RM_SpeciesConcentration2Module(id,...) and RM_GetSpeciesConcentrations(id,...).
As my representative volume is really small in each cell (in the order of 10^-5 litr) I try to scale up the representative volume (to unity) and reactive surface proportional to each other in a way to have representative volume of 1 and get the same results as the small system (if we multiply both solution volume and reactive surface area by a factor we should get the same concentrations in a kinetic reaction), in case of equilibrium it does not make a difference since the system is not time dependent and the amount of reactive component is excessive to reach equilibrium.
previously, I have defined rate expression as below

     RATES 
     anorth
     10 area = PARM(1) # cm^2
     20 k = 10^-7.32*ACT("H+")^1.5 +10^-15.6*ACT("H2O")+10^-17.5*ACT("OH-")^0.33 # mol/s/cm^2
     30 rate = area * k * (1 - SR("Anorthite"))
     40 moles = rate * time
     50 SAVE moles 
     #END 


and updated parm(1) for each reactive cell as the reactive surface area regarding the model specifications (in our model we know exacly the amount of available surface area for each cell).
here the question is that how one can update parm variable in PhreeqcRM?
or how can I defne the amount of reactive surface (or reactive mole) to obtain the desired reactive area?

[dlparkhurst]

I think you will need to write KINETICS_MODIFY blocks for each cell that you want to update. You will need to define all the parms, but you need not change any other kinetic properties for the cell.

When using parallel PhreeqcRM, there is a forward mapping from the transport cells to the chemistry cells. I think you could run all of the KINETICS_MODIFY blocks through all of the workers.

With more work, I you can determine which cell is on which worker (GetEndCell and/or GetStartCell). Then you need to find the IPhreeqc instance for the worker and run just the cells needed for each worker as demonstrated toward the end of the example.

[hamiderfani]

Dear Dr. Parkhurst

I have a question regardin coupling PhreeqcRm with a transport model. In my model I transport all species with my transport model and then run PhreeqcRM for the transport timestep to update species concentrations due to kinetic and equilibrium reactions in each cell, as you know this method is SNIA.
My question is that in order to save time, I can lump some species together and solve transport equations for that pseudo component, as an example i can lum all aluminium bearing species (Al3+, Al(OH)2+, Al(OH)4,...) into [Al]_T which is summation of these ions and solve transport equation for [Al]_T. But my problem is that when I want to update phreeqcRM module concentrations using     "status = RM_SpeciesConcentrations2Module(id, species_c)" I must enter each species concentration (eventhough I can solve transport for each specie instead of lumped species but it is more time consuming) while in solution definition in Phreeqc we enter [Al]_T. is there any way to do this phreeqcRm module instead of defining each species concentration individually?

[dlparkhurst]

PhreeqcRM lets you transport total concentrations of elements, or concentrations of species. There is no way to transport some totals and some species. Unless you are doing electro migration or multicomponent diffusion, you should be able to transport the totals. Otherwise, I think you must transport each species individually.

[hamiderfani]

Dear Dr. Parkhurst
I really appreciate your quick and useful answers. For the time being I have a problem in PhreeqcRM that I really appreciate your help.
recently, I am coupling phreeqcRm with a transport code, which is needed to be run during a long timescale (100 years lets say) as the transport phenomenon is quite slow and diffusion and density deriven.
For this purpose the time steps are quite long, each time step is around 1 month and the phreeqcRm is needed to be run for almost 10,000 gridblocks each time step and we have almost 1000 time steps  totally. I am modeling sandstone rock with some clay minerals and also there is a high possibility for precipitation of carbonate minerals (calcite and dolomite). We also do not want to consider equilibrium condition and we want to use kinetic reactants and their rates. The rates file is attached as well.
My problem is that  the phreeqcRM needs a very very long time in each time step to converge, it gives me the warning for negative values for Al and C, the runtime is very long in a way that makes the code unusable.
I also tried to use -step_divide 2500 ;  -cvode true in kinetics data block but it is still not fast enough and time-step is changed internally by phreeqcRM successively.
It worth mentioning that I am doing the calculation for a representative volume of 1kg of water.
Is there a better way to increase the code speed or to have a robust estimation of the ion concentration and amount of dissolved/precipitated minerals in a reasonable time by PhreeqcRM ?

[dlparkhurst]

Attach your p4w file.

[hamiderfani]

what is p4w file?

I use PhreeqcRm with Fortran, first I initialize the module with desired number of blocks. then I load database and rates file. Then I define the initial condition fluid (based on equilibrium condition with rock minerals), then I do transport and then run PhreeqcRM to update concentrations.
In phreeqcRM I define a kinetic block, which the amount of surface area is defined based on rock composition assumptions and rock properties (to calculate total surface). Kinetics is defined based on reactive surface area of each mineral in gridblock and the rates are also written based on reactive surface area as you can see.

[dlparkhurst]

Sorry, I'm mixing up threads, but I need to see a more complete definition of your reactions, particularly a KINETICS data block that demonstrates the slow calculation.