Using Phreeqc in C++

[GeeqC]

Dear Dr. Parkhurst,

We are trying to link PhreeqC with a C++ code but encountered some problems. In order to run the example for C++ given in the Readme file we followed the instructions in the Readme file (README.txt). First we wanted to try to run the example_advect_cpp, but we failed.

Here is what we have done:
- An Iphreec.sln file was created using Cmake.
- The program was compiled using Visual Studio (ALL_BUILD, RUN_TESTS, and INSTALL) in released mode.
- An *.exe file was created: example_advect_cpp.exe

When the example was run, an error message appeared:
0x000007FEFD63B87D in example_advect_cpp.exe
and
IPhreeqcStop at 0x00000000002CF0E0.

We would be thankful for some support for this procedure.

Best wishes,
GeeqC

[dlparkhurst]

Looks like there is a bug in the INSTALL method or some other problem. We'll have to look at it.

In the mean time, the problem is that the database phreeqc.dat and file ic are not located. You can use ../../../../examples/cpp/advect/ to find the correct files. So modify the source file advect.cpp in the project example_advect_cpp to have

if (this->IPhreeqc_ptr->LoadDatabase("../../../../examples/cpp/advect/phreeqc.dat") != 0) this->EHandler();

and

if (this->IPhreeqc_ptr->RunFile("../../../../examples/cpp/advect/ic") != 0) this->EHandler();

[GeeqC]

Dear Mr. Parkhurst,

Here I would like to follow up on the discussion in the "basics" section on the topic "Converting from molality (mol/kgw) to mg/l for high solubility minerals", where you recommended us to better use the module PhreeqcRM for interaction with a reaction transport model. So, this brings us back to this section here, which is on incorporating Phreeqc in C++ language.

In the meantime, we were able to compile both IPhreeqc and PhreeqcRM for C++. While it was possible to transfer data from the reaction-transport model to Iphreeqc via strings, we do not understand how data transfer works for PhreeqcRM. According to the documentation, it should be possible to transfer entire data arrays which would make the interaction with a reaction-transport model much more efficient. However, we do not know the syntax. In the example given, the data are directly read from a file (advect.pqi) written in Visual Basic. But how can we directly transfer a data array from C++?

Best wishes!
Patrick

[dlparkhurst]

The strategy of using PhreeqcRM is to define a set of solutions and reactants in the "initial iphreeqc" instance. This can be done by either the RunString or RunFile method of PhreeqcRM (in advection_cpp.cpp, the file advect.pqi is read). The chemistry cells for reactive transport are then populated with the method InitialPhreeqc2Module that transfers solutions and reactants to the cells in the "worker iphreeqc" instance where all of the chemistry calculations for the grid are performed. An initial call to RunCells performs any equilibrations in the worker instance and makes the initial conditions from  available to begin transport and reaction. The main time-stepping loop then consists of interacting with the worker instance(s) with GetConcentrations (or GetSpeciesConcentrations, if you are doing multicomponent diffusion), transport calculation, SetConcentrations (or SetSpeciesConcentrations), and RunCells.

There are quite a few methods altogether to control setting units, defining the grid, setting the database, getting lists of elements and reactants, extracting results, controlling parallelization, and so on.

The file Tests/advect_cpp.cpp gives a fairly complete example of many of the most frequently used methods of PhreeqcRM. If you are doing multicomponent diffusion, look at species_cpp.cpp. Documentation of all PhreeqcRM methods is distributed in the directory Doxygen/html and online at https://water.usgs.gov/water-resources/software/PHREEQC/documentation/phreeqcrm/index.html. The reaction-transport codes PHAST and VS2DRT use PhreeqcRM for their chemistry calculations; the source code may or may not be useful examples for you to read.

To answer your specific question SetConcentrations or SetSpeciesConcentrations and GetConcentrations and GetSpeciesConcentrations tranfer arrays of concentrations to and from the reaction module during time stepping. Look for these methods in advection_cpp.cpp or species_cpp.cpp.

[GeeqC]

... it is still not clear to me, how the data can be transferred. Do we have to build the input string in the same way as for Iphreeqc, and including a selected output? How are the concentration vectors read by SetConcentrations transferred to the input string. And how are the selected outputs sent back to the reaction-transport model?

Also, do we have to use RunString or RunCells in this case?

It is furthermore not clear to me, why we need initial conditions. The reaction-transport model already has initial conditions and provides the solutions as input for Phreeqc. And also, how is the "initial Iphreeqc" instance linked to the "worker Iphreeqc" instance?

Thanks for your patience, answering our questions!

[dlparkhurst]

Did you look at advection_cpp.cpp? It demonstrates the sequence of using PhreeqcRM methods.

I will try one more time. Let's say you have 100 cells in your reaction transport model. You need to define the chemistry that is in those cells. Now the chemistry in PHREEQC includes the total moles of H, O, charge, and each element in those cells. Do you know the total amount of H in a solution? No. It takes a PHREEQC calculation to know these quantities.

So, there is a RunFile statement in advection_cpp.cpp. This defines a solution in the "initial" instance of PhreeqcRM. After RunFile, the initial instance knows the complete composition of a solution, including the total moles of H, O, etc. You also define other reactants in the initial instance. Normally, RunString and RunFile are not used again. Then you need to put that solution and reactants (EQUILIBRIUM_PHASES, EXCHANGE, etc) in all 100 cells for which chemistry is going to be run with InitialPhreeqc2Module. Now RunCells will equilibrate the solutions with the other reactants in the cells and the initial conditions are set in PhreeqcRM for all 100 cells. But how will you do the first transport step?

You need to pull the concentrations out of PhreeqcRM to do the transport step. You use GetConcentrations to get an array of concentrations that you are going to transport. That's how you get the total concentration of H in each cell to be able to transport H. Same for all other elements.

You transport, then fill an array with all of the transported concentrations, send them back to PhreeqcRM with SetConcentrations, which updates the concentrations in each cell, and run chemistry with RunCells (not RunString), which takes the new concentrations in each cell and reacts with the remaining reactants each cell. You can then repeat the sequence for the next transport and reaction step.

You can transport and run reactions without using SELECTED_OUTPUT. SELECTED_OUTPUT is used to extract results that you want to examine. If you are looking at a contaminant plume, you could define only the concentration of the contaminant in a SELECTED_OUTPUT definition, and then use the PhreeqcRM methods to retrieve and write the data to file, and then plot it after the simulation is completed.

You seem to be having trouble understanding that PhreeqcRM does all of the work you were doing with strings with IPhreeqc. Everything you were doing with IPhreeqc is done internally (although not necessarily with strings) in PhreeqcRM. You interact only through the 100 or so methods of PhreeqcRM.

[GeeqC]

Thanks for the detailed explanation. While we are still working on the problem, we are wondering whether it is possible to change the path to the database. Where is the path actually defined?

[dlparkhurst]

The following is used in advection_cpp.cpp. It assumes the file phreeqc.dat is in the same directory as the Visual Studio solution, which is the way the tests of PhreeqcRM are run (in the Tests directory).

Code: [Select]

// Load database
status = phreeqc_rm.LoadDatabase("phreeqc.dat");

Alternatively, you can simply give the entire pathname for the database in the call to LoadDatabase.

Note that when you run INSTALL from the PhreeqcRM Visual Studio solution, the databases (and include files, library, Fortran interface code, and other files) are copied to the directory you have chosen when you run CMake (see doc/README.TXT). You could use that directory in your pathname for the database file.

[GeeqC]

Dear David Parkhurst,

In the meantime we got PhreeqcRM to work, and a few simple test runs were successful. But now we encountered a fundamental problem:

The output of the reaction model (and input to PhreeqcRM) consists of the major ion concentrations, total alkalinity, and total C(4). The alkalinity can be well determined based on the conservative ions, according to Wolf-Gladrow et al. (2007) Marine Chemistry 106, 287-300. With these parameters given, the solution speciation should be entirely defined.

However, the set.concentration/get.concentration loop in PhreeqcRM does not include alkalinity but only element concentrations. Running the model with alkalinity defined with the initial conditions did not produce a meaningful outcome. We also tried to fix the pH instead of alkalinity, but also this was not successful.

How can we proceed to correctly define the carbonate system?

Best wishes,
Patrick

[dlparkhurst]

This seems to be a PHREEQC question rather than PhreeqcRM.

You need to determine a SOLUTION definition that correctly defines your solution. For SOLUTION, you need to define two of the three quantities--pH, alkalinity, and total inorganic carbon. The third quantity will be calculated from the other two (particularly, pH will be calculated if you define Alkalinity and C(4)). Once you define SOLUTION, the initial solution calculation will be run and at that point all the moles/concentrations for PhreeqcRM are known, the total H, total O, charge balance, total C, and each of the other elements. If you are happy with the SOLUTION calculation, you are set for PhreeqcRM.

Let's assume this solution is the initial condition for all of your cells. The strategy of PhreeqcRM is to define SOLUTION in the initial IPhreeqc instance of PhreeqcRM; the initial solution calculation will be run, calculating the elemental totals that are needed; InitialPhreeqc2Module will be used to transfer the solution to all the cells in the worker instances; and finally, GetConcentrations will retrieve the element concentrations for use in a transport calculation.

So, it boils down to a PHREEQC SOLUTION calculation. If you get the results that you want from a PHREEQC initial solution calculation, you can use the above strategy to define the solution and transfer it to the workers of PhreeqcRM and then GetConcentrations.

[GeeqC]

Indeed, there was a problem with the definition of the initial solution: I forgot to add "density 1 calc". So this part works fine now.

However, the more fundamental problem is how the alkalinity can be renewed in each round of PhreeqcRM. The reaction-transport model calculates alkalinity based on the diffusion of the conservative ions, but there is no possibility to renew the alkalinity in PhreeqcRM. Instead I would have to renew the oxygen and hydrogen concentrations, but these are not known. They are calculated by PhreeqC.

Alternatively, oxygen and hydrogen concentrations could possibly be calculated as the sum of all species containing these elements. But this would be rather complex and diffusion coefficients for many species, such as ion pairs, may not be known. Moreover, charge balance may also change due to diffusion. I wonder whether there could be a more simple way to solve this problem.

[dlparkhurst]

I do not know whether you are transporting each aqueous species as in the species_f90.F90 or total element concentrations as in advection_f90.F90. It depends on whether diffusion is the main transport mechanism (you would transport species with diffusion coefficients for each species) or advection is the main transport mechanism (transporting total elements is sufficient).

In either case, you would not need to transport alkalinity. You would do transport either with individual species (including water) or with total elements (including H, O, and charge), send the concentrations back to PhreeqcRM, and the next RunCells calculation would calculate the new alkalinity after reaction. If there are no other reactions, RunCells would calculate the alkalinity due only to the effects of transport. You can obtain the alkalinity for the cells by processing results from SELECTED_OUTPUT (-alk) or USER_PUNCH (PUNCH ALK).

If you know the total concentrations of species in a solution (and the charge), then you can calculate the total concentrations of elements in solution. If you know the total concentrations of elements (H, O, charge, C, others), that is sufficient for PHREEQC to calculate pH, TDIC, and alkalinity.

[GeeqC]

Indeed, we are simulating diffusion, for which the transport of different species matters. However, the main array "c" in PhreeqcRM only contains totals, not species. How can we transfer species back to PhreeqcRM?

I assume that all species would have to be summed up to total elements. But there are so many species, including all the ion pairs, that it would be extremely complex to simulate diffusion of all species correctly.

For simplicity I used total elements and assumed that diffusion is dominated by the most abundant species. But then I run into problems when I have to calculate total H, O, and charge. Would there be a better way to proceed?

[dlparkhurst]

You have two choices: you can either (1) transport all of the species using GetSpeciesConcentrations and SpeciesConcentrations2Module (and other methods) with diffusion coefficients for each aqueous species as demonstrated in the example species_cpp.cpp, or (2) you can transport total concentrations of each element (GetConcentrations and SetConcentrations) using a single diffusion coefficient as demonstrated in the example advection_cpp.cpp.

Sorry if it is hard. The pitzer.dat database has very few aqueous species, but also contains few elements (and its numerical method is slower). phreeqc.dat has diffusion coefficients for many, but not all species.

[GeeqC]

It makes a lot of sense to use species for diffusion calculation. And probably it may not be necessary to transport all the minor species. For example, diffusion of H+ is probably negligible. The correct H+ concentration will be calculated in each round with the speciation.

I am not sure, however, if diffusion of H2O matters, and how we should set the boundary conditions for H2O. We are already calculating the H2O content using "density 1 calc".

[dlparkhurst]

You can decide how you want to implement the diffusive transport calculation, but if you are developing a general-purpose program, there would be cases where H+ diffusion would be critical.

Water should diffuse according to its activity.

[GeeqC]

Is it possible to calculate the species concentrations (not activities)? Because to my understanding the diffusion coefficients available (e.g. in Boudreau, 1997) are for concentrations under specific P, T, and salinity conditions.

[dlparkhurst]

You can look at the documentation for the TRANSPORT data block in the PHREEQC manual to see the equations used for multicomponent diffusion. I think you need to be careful of shortcuts and simplifications.

No PhreeqcRM method returns activity. PhreeqcRM GetSpeciesConcentrations returns concentrations of species in mol/L, and GetSpeciesLog10Gamma returns the activity coefficients. You need to be  careful because the activity coefficients are for molality rather than molarity.

[GeeqC]

Yes, the concentrations in mol/L are fine. We calculate the diffusion coefficients based on Boudreau (1997).
But what I do not understand is that you mentioned water activity. Why is it "activity" for H2O, but "concentration" for all other species?

[dlparkhurst]

I think it is activity for all species. Again, check the equations and references in Phreeqc Version 3 manual on TRANSPORT.

[GeeqC]

Dear Mr. Parkhurst,

It works! Thank you so much for guiding us through this procedure.

While playing around with the parameters, we encountered one problem:
We are simulating seawater in the deep ocean at 5000 metres waterdepth. If we use a pressure of 500 atm, the sum of all sulphur species and ion pairs becomes larger than the total element contentration. The difference is on the order of 0.5 mmol/L. This only occurs in PhreeqcRM, while the correct concentrations appear in Phreeqc with User Interface. What could be the cause of this discrepancy?

Best wishes,
Patrick

[dlparkhurst]

My guess is that you are looking at different units, perhaps mol/kgw vs mol/L.

If you print selected output values in PhreeqcRM, The sum of species, using the function MOL (mol/kgw) will equal the total concentration, using TOT (mol/kgw).

Concentrations are converted for the transport code to mg/L, mol/L, or mass fraction, as specified by the PhreeqcRM method SetUnitsSolution. Results from PHREEQC are generally in mol/kgw, although you can use Basic functions GFW (g/mol), RHO (solution density, kg/L) TOT("water") (kg water), and SOLN_VOL (solution volume in liters) to print selected output values in mg/L, mol/L, or mass fraction. The following calculates the distribution of sulfur species in mol/L and compares the sum of species to the total sulfur in solution. Note, in this example, the stoichiometric coefficient of sulfur in the aqueous species needs to be taken into account because of the species Fe(SO4)2-.

Code: [Select]

SOLUTION
pH   7 charge
S(6) 1
Fe(3) 1
USER_PRINT
10 t = SYS("S", count, name$, type$, moles)
20 PRINT "Species                  mol/L           S coef"
30 FOR i = 1 TO count
40    if (type$(i) <> "aq") THEN GOTO 110
50    n$ = SPECIES_FORMULA(name$(i), count_e, elt$, coef)
60    FOR j = 1 TO count_e
70      if (elt$(j) <> "S") then goto 100
80      sum = sum + MOL(name$(i)) * coef(j)
90      PRINT PAD(name$(i),17), STR_E$(MOL(name$(i))*TOT("water")/SOLN_VOL, 15, 4), coef(j)
100   NEXT j
110 NEXT i
120 PRINT "Sum species, mol/L: ", sum*TOT("water")/SOLN_VOL
130 PRINT 'Total S, mol/L:     ', TOT("S")*TOT("water")/SOLN_VOL
END


[GeeqC]

... In the meantime we did further tests to locate the problem. Here some specifics:

1.) SetUnitsSolution was set to mol/L

2.) The input string was set as follows:

-units mmol/l
-density 1029 calc
-S(6)  28
SELECTED_OUTPUT
-totals S(6)

The string (RunString) was still calculated correctly.

3.) After RunCells the "totals S(6)" in the selected output were still correct (28 mmol/L), even at high pressure.

4.) The sum of the main S(6)-species (SO4-2 + Ca(SO4) + Mg(SO4)) read from the "c_species" array was still correct (i.e. 28 mmol/L). Only at high pressure (500 atm) the sum of species became ca. 28.5 mmol/L (i.e. larger than the total element concentration).

In conclusion: the units seem to be correct, as otherwise we should already see wrong results at low pressure. The question is what the pressure has to do with it.

Best,
Patrick

[dlparkhurst]

Solution volume is not constant with pressure, so concentration in mmol/L (or mg/L) is not constant with pressure. Molality would be more constant, although speciation could change with pressure, so, the mass of water could vary slightly if hydrolysis reactions were pressure dependent. Mass fraction of S should be constant with pressure. In all cases, density is not constant with pressure, and your transport model would need to take that into account if pressure varies significantly.

In the attached file, the solution volume for a 28 mmol/L solution scaled to 1 kilogram of water is 1.00334 L at 1 atm pressure. However, the solution is calculated to be 28.6 mmol/L with a volume of 0.98201 L at 500 atm. Moles of sulfate is the same in both solutions, but volume has changed. Molality is 28.19 mol/kgw in both solutions

Code: [Select]

SOLUTION
-units mmol/L
-density 1.029 calc
pH 7 charge
Na 56
S(6) 28
USER_PRINT
10 PRINT "SO4, mmol/L: ", 1000*TOTMOL("S(6)")/SOLN_VOL
END
USE solution 1
REACTION_PRESSURE
500
END


[GeeqC]

Yes, this makes sense. Thereby, the boundary conditions also need to be adjusted to the high pressure.

[GeeqC]

PhreeqcRM runs now very well. However, I noticed that some reduced species, such as NH4+, HS-, and CH4 become depleted at each equilibration step. This occurs also when the transport function is switched off. I assume it has something to do with the equilibration. Is it possible that these species undergo oxidation? The only oxidant in the solution would be SO42-.

Is it possible to turn the redox reactions off?

[dlparkhurst]

You could turn redox reactions off using the method of example 9, where Fe is split into two, non-interacting elements.

However, I doubt you want to do that. You should figure out why the redox reactions are occurring. If N(-3) is oxidized, then something must be reduced, and from your statement, it is not S(6) because S(-2) is also being oxidized. You may have an oxidant that you are not considering (O(0), Fe(3), or N(5), for example), or you may have defined a REACTION or KINETIC formula incorrectly. As an example of a formula that would oxidize, NaO2 would be equivalent to NaO0.5 + 0.5O2 (or essentially NaOH + 0.5O2), although an erroneous formula could be much more subtle.

If you can extract a calculation as a PHREEQC input file (perhaps using DUMP before one of the calculation), I will explain the redox reactions.

[GeeqC]

This is the input file:

SOLUTION 1
   -temp   10
   -units   mmol/L
   -density   1   calc
   Alkalinity   2.3
   C(4)      2.1
   C(-4)   1
   N(-3)   1
   Ca      10
   Mg      50
   Na      470
   K      10
   Cl      550
   S(6)      28
   S(-2)   1
   Si      1
   -water 1 # kg

There are no further electron acceptors used in this simulation. I do not use kinetics, just equilibration. In this case, should the total concentrations of N(-3), C(-4), and S(-2) not remain constant?

Then there should be another reason why NH4+, CH4, and HS- are disappearing. What could it be?

[dlparkhurst]

SOLUTION allows redox disequilibrium.  When you define a valence state like S(6), that is the total sulfate concentration, regardless of redox equilibrium. When you define a reaction, here MIX, then all of the redox states adjust to equilibrium.

One reaction that happens to produce redox equilibrium is the production of N2(aq). It is often reasonable to assume N2(aq) does not react, in which case, you can disable the reaction by redefining N2(aq) with a smaller log K, as I have done in this example.

Code: [Select]

SOLUTION_SPECIES
2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
#-log_k 207.08
      -log_k  0
-delta_h -312.130 kcal
-dw 1.96e-9
-Vm 7 # Pray et al., 1952, IEC 44. 1146

SOLUTION 1
   -temp   10
   -units   mmol/L
   -density   1   calc
   Alkalinity   2.3
   C(4)      2.1
   C(-4)   1
   N(-3)   1
   Ca      10
   Mg      50
   Na      470
   K      10
   Cl      550
   S(6)      28
   S(-2)   1
   Si      1
   -water 1 # kg
END
MIX 1
1 1
END
[/code

Another reaction is the oxidation of C(-4) by S(6), where C(-4) basically disappears, S(6) decreases, and S(-2) increases. I don't know if you want to disable this reaction of not; you could argue either way. If you want C(-4) not to react, you can define the methane as Mtg instead of C(-4). Mtg is not defined as part of the C (carbon) aqueous species. Mtg has only one aqueous specie and cannot be oxidized, although it could increase or decrease by REACTION or KINETICS. Note that CH4(aq) can still form if, for instance, organic matter (CH2O) reacts with your solution; Mtg will not change but C(4) could be reduced to CH4(aq).

[GeeqC]

Sorry, I meant to say that equilibration of redox states should not occur in this simulation. Only speciation.

N(0) is not defined in this system, thus, N2 should not form. The total concentrations of N(-3), C(-4), and S(-2) should remain constant.

The question is, what could cause the disappearance of NH4+, CH4, and HS-?

[dlparkhurst]

I would need to have your database, or you can run the script of my previous post and attach the output.

[GeeqC]

Sorry again. I did not use the correct Phreeqc-language. What I meant to say is that N(0) is not given in the input file. The database is sit.dat.

Redox reactions should not occur for N-species, because only one redox state occurs in the input file.
If I run Solution 1 as above (without kinetics) in the regular PhreeqC user interface, then all species are there. The only species appearing are NH4+, NH3, and some ion pairs. I considered NH4+, NH3, and Mg(NH3)+2 in the simulation. Why are they being consumed when itererating through the PhreeqcRM module?

[dlparkhurst]

Sit.dat has both N(-3) and N(5). Assuming that transport is not responsible for a decrease in N(-3), then it must be reacting through exchange, surface, kinetics, etc., or it is oxidized to N(5). The reaction calculations of PHREEQC should conserve mass. So, is N(5) produced?

I suggest that you remove any reactions like equilibrium_phases, kinetics, exchange, and surface, and then look at N(-3) and N(5) concentrations.

There are a lot of moving parts, which makes debugging a reactive transport code difficult. If you can isolate a PHREEQC script that demonstrates the reactions, I will look at it.

Usually, you run an equilibration step before transport time stepping. Does this reaction occur during the initial equilibration or during transport?

You can also define your initial solutions to be in redox equilibrium (through an initial PHREEQC calculation) and constant over the domain (and boundary conditions) with no additional reactions. Reaction following transport should produce no changes.

[GeeqC]

Thanks for your patience! Now I could find the problem: plotting total N(5) as selected output, indeed, it builds up while N(3) is depleted.

So in conclusion: PhreeqcRM is different from PhreeqC in that it changes redox states. In my case, I would like to turn off all redox reactions. This is because under Earth surface conditions most of these reactions are inhibited, and they only occur via enzymatically controlled kinetics. These are calculated externally, using metabolic rate laws.

Is it necessary to turn off each redox reaction individually (as in your earlier example or in Example 9), or is there a way to turn off redox reactions generally?

[dlparkhurst]

You need to figure out why N(-3) is oxidized. I think you must have O2 in the system, an oxidizing reactant, or an error in your transport. PHREEQC and PhreeqcRM are identical when given identical inputs.

Note that with PhreeqcRM you must transport H, O, charge, and the totals of each element. Alternatively, you can transport H2O, excess H (H - 2H2O), excess O (O - H2O), and charge. The relatively small amounts of oxygen and hydrogen not in water are critical to the calculation of redox states.

[GeeqC]

In this case species are transported, not elements (because different species have different diffusion coefficients). I have now included H+, OH-, and H2O, but this did not change the outcome.

The transport is calculated externally. The transport model works fine if PhreeqcRM is switched off. Also PhreeqcRM works fine if the transport model is switched off. Possibly, I am missing some essential species.

[dlparkhurst]

You must transport all species. The method GetSpeciesNames provides a list of all aqueous species that you need to transport. There are prerequisite method calls, so look at the documentation. GetSpeciesZ returns the charge on the species; an incorrect charge or a missing species could cause your problem.

You should use charge balanced solutions as well, which is probably a requirement for multicomponent diffusion.

The example species.cpp in the PhreeqcRM distribution demonstrates the use of species for transport, although it does not calculate multicomponent diffusion.

PHREEQC can calculate multicomponent diffusion, so that could be a way to check your scheme.

[GeeqC]

How can I charge-balance the solution? The input file only has totals, the charges are calculated by Phreeqc.

[dlparkhurst]

Make sure you charge balance all of your solutions in the initial PHREEQC run.

Usually, the only reaction that can introduce charge imbalance into a solution is SURFACE when not using -donnan. So, if you are using SURFACEs, make sure you use -donnan (or -edl).

[GeeqC]

In the meantime, I essentially did everything as suggested:

- The external source and sinks were charge balanced with H+ (respectively, they were entirely switched off).
- The diffusion equations were added for all species that occur in relevant amounts, including H+, OH-, and the ion pairs.
- As upper boundary condition I used the species concentrations from the initial input file. I had to calculate them beforehand using the Phreeqc User interface. However, the precise species concentrations slightly deviate in the initial solution calculated by PhreeqcRM.
- As the lower boundary condition I assumed a zero gradient (Dirichlet-Boundary condition) for each species
- The charge balance: this is still a problem, because phreeqc is calculating the charge balance based on speciation. I can vary the totals in the input file and try to minimize the charge by trial and error. However, a precisely charge-balanced state cannot be reached in this way. Is there an alternative way, so that Phreeqc charge-balances the solution automatically?
- I ran the species output to check all the species concentrations. Everything looks fine, i.e. no more buildup of O2 or NO3-.

The outcome is that all species concentrations go offscale. Something is wrong. I tested the diffusive transport alone, but this part works fine. The problem occurs when I turn on PhreeqcRM.

[dlparkhurst]

You can use the "charge" option in SOLUTION with PHREEQC to produce a perfectly charge balanced solution.

You should closely follow the sequence of method calls in the example species_cpp.cpp (or species_f90.F90) in the Tests directory of the PhreeqcRM distribution. Generate your initial and boundary conditions with an initial run of PHREEQC (where you charge balance all the solutions) followed by InitialPhreeqc2Module to populate your initial cell concentrations. Use InitialPhreeqc2SpeciesConcentrations to get boundary condition concentrations. I would transport every species, whether you think it is important or not.
 
After setting up initial and boundary conditions, the first check would be to use SpecieisConcentrations2Module followed by GetSpeciesConcentrations, with no call to RunCells to see if you get back the species concentrations that you put in. Then, without using any reactants (EQUILIBRIUM_PHASES, etc) in the cells, add the call to RunCells. If you do not get back your initial concentrations, you probably have a problem with units conversion.

[GeeqC]

As you suggested, I did the following modifications:

1.) "charge" was added to Cl (because chloride is very abundant and it does not matter if its concentration is slightly adapted to balance the charge).

2.) "InitialPhreeqc2Module(ic1)" is used to populate the array with the initial concentrations.

3.) I did not use "InitialPhreeq2SpeciesConcentrations" because boundary conditions are not needed in Phreeqc. I am only running a speciation, while the transport is calculated externally.

4.) All species are now transported.

Overall, this improved a lot the simulated profiles, but there still are a few problems. Even if "charge" is assigned to Cl in the initial solution, the charge balance changes over time. If the pressure changes with depth, I still get strange values. The charges are no longer balanced throughout the profile. How can I balance the charges at all depths if pressure and temperature are changing?

[dlparkhurst]

If you have boundary condition concentrations in your transport, I think you should use PHREEQC (InitialPhreeq2SpeciesConcentrations) to define these concentrations. You need to have exact (charge balanced) concentrations of all the species at the boundary.

Charge balance for a solution in a cell of PhreeqcRM is determined by the charge on the species times the charge of the species summed over all species, as calculated from the array given by SpeciesConcentrations2Module. So, if charge balance is changing, it is because the transported species no longer produce charge balance. The charge balance calculated in this way will not change during any reaction calculation occurring in RUN_CELLS.

[GeeqC]

Then I do not understand why to use InitialPhreeq2SpeciesConcentrations. The boundary conditions for the transport calculation are the same as the concentrations calculated for the initial solution and can be extracted using GetSpeciesConcentrations.

[dlparkhurst]

I just wanted to make sure that you are using PHREEQC-derived concentrations for the boundaries. Often the boundaries are different from the initial conditions, so InitialPhreeqc2Concentrations would be the way to get exact, charge-balanced concentrations for these non-initial-condition solutions.

[GeeqC]

I see. This could be useful to set fixed lower boundary conditions that are different from the initial solution (so far I used zero flux).

But how would I do that? I can see how InitialPhreeq2SpeciesConcentrations is used in the example species.cpp. But would it not be necessary to define a new initial solution for that lower boundary?

[dlparkhurst]

There are three sets of IPhreeqc instances in PhreeqcRM. A set of workers, an initial-phreeqc instance, and a utility instance. The set of cells corresponding to the transport model are stored in the workers. The idea is to run a PHREEQC file in the initial-phreeqc instance, and then extract the definitions that you want as initial conditions and boundary conditions. Initial conditions get transferred to the workers with  SpeciesConcentrations2Module, and boundary conditions are retrieved with InitialPhreeq2SpeciesConcentrations and processed by the transport model.

So, you can define multiple SOLUTIONs in the PHREEQC input file, some of which are used for initial conditions and some used as boundary conditions. Simplest case, SOLUTION 1 defines initial conditions, SOLUTION 2 is a boundary condition. SpeciesConcentrations2Module sets SOLUTION 1 as the initial conditions, InitialPhreeq2SpeciesConcentrations gets concentrations of species in SOLUTION 2 for boundary conditions.

[GeeqC]

In the meantime, I tried to use InitialPhreeqc2SpeciesConcentrations. It works, but I do not see the advantage over GetSpeciesConcentration. Both commands can be used to read out the species concentrations from initial solutions, both Solution 1 and Solution 2.

This exercise revealed further issues:
- If I have an array with 20 depth intervals, is it possible to use only one initial solution? The charge balance is then the same for all of them, even if the temperature and pressure changes with depth?

- If I use Solution 1 as initial condition and Solution 2 as lower boundary condition, they have different charge balances. Is it not a problem that at the lowest depth interval, the solution jumps from one charge balance to the other?

[dlparkhurst]

>In the meantime, I tried to use InitialPhreeqc2SpeciesConcentrations. It works, but I do not see the advantage over GetSpeciesConcentration. Both commands can be used to read out the species concentrations from initial solutions, both Solution 1 and Solution 2.

I will try again. The "Initial Phreeqc" IPhreeqc instance is intended as work space to define solutions for initial and boundary conditions. The "worker" Iphreeqc instance(s) are intended to contain the data for the cells during transport simulations. So, the plan is to define some solutions in the Initial Phreeqc instance and use them to define the initial conditions for the workers (SpeciesConcentrations2Module). You can also use InitialPhreeqc2SpeciesConcentrations to get data for boundary conditions for your transport code. If boundary conditions don't change, then you are done with the Initial Phreeqc instance.

GetSpeciesConcentrations gets concentrations from the worker instances (InitialPhreeqc2SpeciesConcentrations gets concentrations from the Initial Phreeqc instance). So, once you have populated the cells in the workers (SpeciesConcentrations2Module), you would use GetSpeciesConcentrations to get the initial conditions to start your transport simulations. After that for transport steps, you should use SpeciesConcentrations2Module to transfer the transported concentrations to the module (workers), run reactions in the cells (RunCells), and retrieve reacted concentrations from the workers (GetSpeciesConcentrations).

It may or may not make a difference, but SOLUTION 1 is just a solution in the Initial Phreeqc instance, but SOLUTION 1 in a worker is associated with cell 1 of the reactive transport model. (Note that with parallelization, there are multiple workers; SOLUTION 1 could exist in all of the workers, but is only meaningful in the worker that is assigned cell 1.)


>This exercise revealed further issues:
>- If I have an array with 20 depth intervals, is it possible to use only one initial solution? The charge balance is then the same for all of them, even if the temperature and pressure changes with depth?

It is possible to assign one solution from the Initial Phreeqc instance to all of the cells in the workers. The charge balance would be the same for all cells. For PHREEQC, any change in temperature or pressure will not change the charge balance.

>- If I use Solution 1 as initial condition and Solution 2 as lower boundary condition, they have different charge balances. Is it not a problem that at the lowest depth interval, the solution jumps from one charge balance to the other?

I would charge balance both solutions in the Initial Phreeqc instance using the "charge" option within the SOLUTION definition. SpeciesConcentrations2Module will then populate the cells with charge balanced solutions. Same for the boundary condition.

If you do this, the only source for charge imbalance is your transport module; PHREEQC maintains the charge balance of the solution in all reactions, including temperature and pressure change (except possibly SURFACE). When using species transport, the charge balance for a cell will be defined by the sum of z*m(i) of the transported species. See TRANSPORT in the Phreeqc manual for how charge balance is maintained in multi-component diffusion.

[GeeqC]

Thanks for the detailed explanation.

In summary: InitialPhreeqc2SpeciesConcentrations gets the concentrations from the initial solution, whereas GetSpeciesConcentrations gets the concentrations from the array of workers.

This means that if I define SOLUTION 2 as boundary condition for the lower domain boundary at a specific Depth z, I would have to calculate this solution for the correct pressure and temperature at this depth. Otherwise the species concentrations of the boundary condition would be wrong?

Alternatively, I could use constant flux as boundary condition. However, I do not know the fluxes of each single species at the lower domain boundary. Only the fluxes of the total elements are known. Would this still allow me to use constant flux as lower BC?

[dlparkhurst]

Yes, you would need to use the correct temperature and pressure for the boundary condition solution. If temperature and pressure are changing during the simulation, you can always go back to the initial phreeqc instance and use SOLUTION 2 with REACTION_TEMPERATURE and REACTION_PRESSURE and then retrieve the adjusted species concentrations.

Not sure what to tell you about how to convert total fluxes to species fluxes. It actually does not matter to PhreeqcRM. When you use SpeciesConcentrations2Module, PhreeqcRM simply sums up the totals of elements from the species for a cell and redistributes the species in the RunCells calculation, so as long as the species sum to the same totals, the result would be the same.

But it does matter for the transport simulation if an element is more in a charged aqueous species or an uncharged species. (I assume this is a diffusive system and you really need to transport species rather than totals.)

So, I do not have the right answer. Here are some thoughts. In Phreeqc, you could simulate a flux with a KINETIC reaction with the formula containing the elements; but I don't see how it would work to transfer the flux to the Phreeqc side rather than the transport side. You could use the initial Phreeqc instance (or perhaps the utility instance), where you take the composition of the boundary cell, add a REACTION for some period of flux, and extract the changes in species concentrations for that period of time. Perhaps that would give you an appropriate distribution of species for the flux.

[GeeqC]

I realize that there are still some issues related to charge balance. The calculation of diffusion based on chemical potential, as described in the manual, would result in a perfect charge balance. However, I have been using Fick's model so far, in which case a charge disbalance occurs.

Surprisingly, however, the simulation runs just fine, since I have charge balanced the initial solution. There are no longer strange artefacts with oxidized species, despite the fact that there is a charge disbalance of 1 meq/l or more. This makes me wonder: what is the difference between charge balance of the initial solution and charge balance of the solution array in the workers?

In other words: how is it possible that the model is "immune" to charge disbalance in the workers? And besides, it is, apparently, possible that electric currents are running through the sediment (e.g. Risgaard-Petersen et al., 2012). Should not the electrons also be considered as a species that is transported?

[dlparkhurst]

I don't know what happened when you used charge imbalanced boundary solutions, but I suspect the following illustrates the difference. In one case the charge imbalance of the solution is not adjusted for the change in concentrations, and in the other, the charge imbalance is adjusted. The first produces a  presumably erroneous pH of 3, and the second produces a pH of 7.

When you transport with Fick's law, you end up with a set of species concentrations for a cell that do not charge balance. However, PhreeqcRM calculates the charge imbalance and updates the charge imbalance of the cell solution as in solution 2 of the example. It is essentially generating a ficticious, inert ion that exactly accounts for the charge imbalance.

I suppose this more or less works for your simulation, although the concentrations definitely contain errors and you may get criticized by reviewers. I would certainly run the same simulations assuming every species has the same diffusion coefficient (either by setting equal diffusion coefficients for species, or transporting totals instead of species) for comparison. Better yet would be to implement true multicomponent diffusion with exact charge balance.

Code: [Select]

SOLUTION 1-2
-units mmol/kgw
pH 7 charge
Na 1
Cl 1
END
SOLUTION_MODIFY 1
-totals
Cl 1.5e-3
Na 1e-3
END
SOLUTION_MODIFY 2
-cb  -0.0005
-totals
Cl 1.5e-3
Na 1e-3
END
RUN_CELLS
-cells 1-2
END


[GeeqC]

While we are trying to implement charge balanced transport, there is also a problem with aluminium:

The Al species go to zero while all other species are simulated correctly. I checked several times for errors in the transport code, but everything seems to be fine. Could there be an other reason for the problems with Al?

[dlparkhurst]

Solubility of Al minerals can be quite small, especially at neutral pH. PhreeqcRM does censor an element total at some concentration (somewhere between 1e-13 and 1e-25, I would have to check). So, check whether the total concentration is in this range.

[GeeqC]

The concentrations of Al are around 1E-5 mmol/l in the first run, but then they drop to zero (except near the top where the concentration is fixed by the boundary condition).

The aluminium is needed in order to calculate the saturation state of clay minerals.

[dlparkhurst]

Then you need to figure out why Al is going to zero.

Skip RunCells and see if everything is OK. I don't see why PhreeqcRM would be a problem at those concentrations.

[GeeqC]

If RunCells is swiched off, everything runs fine.

[dlparkhurst]

Turn off all the reactions (EQUILIBRIUM_PHASES, EXCHANGE, etc), and check the error code from RunCells.

[GeeqC]

All reactions are turned off. There is no error message. The simulation runs all the way to the end, just the Al concentrations go to zero.

[dlparkhurst]

Does Al go to zero after one reaction step?

Sorry, I don't have an idea of what is happening, unless you are missing an aqueous species. I assume you are getting the list of species from PhreeqcRM.

I think I would try a different database first, just because that should be easy. Then I would turn on and print out the PHREEQC output for each cell. That should tell you what PHREEQC is calculating.

[GeeqC]

It took a bit longer, but several errors were found. The Al and Fe are now calculated correctly. All transport is now charge balanced.

Now we can come back to the boundary conditions: We define SOLUTION 1 as initial condition. The upper boundary condition is equal to SOLUTION 1. For the lower boundary condition we define a SOLUTION 2.

Next, the array c_species has to be populated correctly. According to the example in "species.cpp" this is done by: status = phreeqc_rm.InitialPhreeqc2Module(ic1);
But how does the program know that SOLUTION 1 should be used?

The same is also unclear for the boundary condition. How does phreeqc know that SOLUTION 2 should be used for the boundary condition in the command:
status = phreeqc_rm.InitialPhreeqc2SpeciesConcentrations(... );

[dlparkhurst]

ic1 contains index numbers for SOLUTIONS, EQUILIBRIUM_PHASES, etc. You need to specify 1 in the solutions slot to use SOLUTION 1 for a cell. Here is the code from the advection_cpp.cpp example that sets the initial conditions.

Code: [Select]

for (int i = 0; i < nxyz; i++)
{
ic1[i] = 1;              // Solution 1
ic1[nxyz + i] = -1;      // Equilibrium phases none
ic1[2*nxyz + i] = 1;     // Exchange 1
ic1[3*nxyz + i] = -1;    // Surface none
ic1[4*nxyz + i] = -1;    // Gas phase none
ic1[5*nxyz + i] = -1;    // Solid solutions none
ic1[6*nxyz + i] = -1;    // Kinetics none
}

It is up to you to specify the boundary conditions for the your transport calculation. You can get the concentrations with InitialPhreeqc2SpeciesConcentrations(... ), but you will need to install them as boundary conditions.

[GeeqC]

How can I get the species concentrations for Solution 2?

[dlparkhurst]

Code: [Select]

std::vector<double> bc_conc;
std::vector<int> bc1;
bc1.push_back(2);                      // solution 0 from Initial IPhreeqc instance
status = phreeqc_rm.InitialPhreeqc2SpeciesConcentrations(bc_conc, bc1);

The concentrations will be in the array bc_conc. To get the species names and order, you need to do the following soon after you run the input file:

Code: [Select]

status = phreeqc_rm.SetSpeciesSaveOn(true);
int ncomps = phreeqc_rm.FindComponents();
int npecies = phreeqc_rm.GetSpeciesCount();
const std::vector<std::string> &species = phreeqc_rm.GetSpeciesNames();


[GeeqC]

In the line "status = phreeqc_rm.InitialPhreeqc2SpeciesConcentrations(bc_conc, bc1);" an error message shows up: "No matching member function for call to 'InitialPhreeqc2SpeciesConcentrations'"

Where is it indicated that bc_conc refers to SOLUTION 2 (instead of SOLUTION 1)?

[dlparkhurst]

Come on. Did you look at the documentation for InitialPhreeqc2SpeciesConcentrations?

bc1.push_back(2);

This line puts a 2 in the vector bc1; this 2 refers to solution 2. InitialPhreeqc2SpeciesConcentrations(bc_conc, bc1) means put concentrations in bc_conc for each solution number found in bc1. In this case, the nspecies items of vector bc_conc are the species concentrations from solution 2 in the initial phreeqc instance.

I don't know why you got an error message. Either you misspelled the method or perhaps have the wrong data types for the vectors. InitialPhreeqc2SpeciesConcentrations is used in the test case species_cpp.cpp, so if the test cases for PhreeqcRM run, then the method is working correctly.

[GeeqC]

Thanks for your patience. Your explanations were very helpful. The boundary conditions are now set correctly.

As a further element, we would like to include mineral reactions. In our case, the reactions are far from equilibrium and so I assume we would have to use the kinetics data block. However, unlike in batch reactions, where the kinetics is integrated over a time interval, this would not be needed for the reaction transport model. The model is already time transient, so we only need to calculate the kinetics at each single time step.

In the user manual it is described how to calculate the time-integrated reactions, but I do not understand how to include kinetics at a single time step in the RT-model.

[dlparkhurst]

A call to RunCells in PhreeqcRM represents the chemical reactions over the time step of the transport model, so you would use KINETICS just like any other situation. KINETICS will be integrated over the time specified by SetTimeStep.

[GeeqC]

I think in our case the integration is not needed, because in the numerical reaction-transport model dt is already the infinitesimal time step:

Rate = (c(t+dt)-c(t))/dt
dc =  (c(t+dt)-c(t))/dt * dt = c(t+dt) - c(t)

Is there a way to switch the integration off?

[dlparkhurst]

You can try using a lower order Runge Kutta method of KINETICS. If your rate is really nearly constant over the interval, then it will do fewer evaluations and run more quickly. However, it the rate is not constant, it will automatically use a higher order, and I am not sure that it will revert back to a lower order.

You might be able to evaluate the rate at each time step in a USER_PUNCH block and apply it in SetConcentrations at the start of the next time step. You need to be careful of units. I wouldn't try it until you are in the last stages of optimizing your program. I think you would be better off to use the parallelization with a lot of processors. At some point, with enough cells and enough processors, the transport calculation takes longer than the chemistry calculation.


Line 9: -runge_kutta ( 1 , 2 , 3 , or 6 )

( 1 , 2 , 3 , or 6 )--This parameter affects integration by the Runge-Kutta solver. It designates the preferred number of time subintervals to use when integrating rates and is related to the order of the integration method. A value of 6 specifies that a 5th order embedded Runge-Kutta method, which requires six intermediate rate evaluations, will be used for all integrations. For values of 1 , 2 , or 3 , the program will try to limit the rate evaluations to this number. If the -tolerance criterion is not satisfied among the evaluations or over the full integration interval, the method will automatically revert to the Runge-Kutta method of order 5. A value of 6 means that the 5th order method will be used exclusively. Values of 1 or 2 are mainly expedient when it is known that the rate is nearly constant in time. Default is 3 .

[GeeqC]

We would like to simulate kinetics of clay mineral dissolution and precipitation. These reactions are essentially limited by aluminium (an partially Fe), which show a very low solubility. The Runge-Kutta method may prevent the reactions from overshooting within a single time step. This effect would slow down the reaction. However, in reality the depletion of Al due to precipitation of one mineral would be compensated by the dissolution of another mineral. Is there a way that phreeqc can resolve these processes within a time step?

Alternatively, the dissolution and precipitation could just be calculated as source and sinks in the reaction-transport model, so that the dissolution and precipitation are stoichiometrically coupled through aluminium.

Also, I am wondering if there is an example for an application of the kinetics function in PhreeqcRM. Because in the example advection.cpp the kinetics function is switched off:

ic1[6*nxyz + i] = -1;    // Kinetics none

and no kinetics data block is given in the input file "advect.pqi".

[dlparkhurst]

If you have mineral sources of Al in EQUILIBRIUM_PHASES, then these aluminum-bearing minerals will dissolve/precipitate in response to a KINETICS reaction. Alternatively, you could have other KINETICS reactions that add or remove Al.

KINETICS work the same way as EQUILIBRIUM_PHASES or other reactions. You would define a kinetic reaction in the advection.pqi file, and distribute it to cells with InitialPhreeqc2Module method. If you define KINETICS 10 in the pqi file that you want in every cell, then put it in the input file and set ic1[6*nxyz + i] = 10. Examples 6, 9, and 15 in the Phreeqc manual demonstrate the use of KINETICS; PhreeqcRM works identically.

Each cell of PhreeqcRM is like a batch reaction in PHREEQC. You should test out your reaction scheme with PHREEQC before adding the complication of transport.

[GeeqC]

Thanks! The examples are helpful, but I will have some further comments to this problematic.

Before I come to this, I would like to address another issue: There is a further problem with the boundary conditions. If I set the temperature and pressure gradient for the array, using:

phreeqc_rm.SetTemperature(temperature);
phreeqc_rm.SetPressure(pressure);

the temperature and pressure are not correctly set for the boundary condition. This results in a kink of the simulated concentration profiles in the lowest depth interval. I tried to fix the problem by defining the pressure and temperature for Solution 2:

oss << "USE solution 2 " << "\n";
oss << "REACTION_PRESSURE 2" << "\n";
oss << pressure[N-1] << "\n";
oss << "REACTION_TEMPERATURE 2" << "\n";
oss << temperature[N-1] << "\n";
oss << "END\n";

But this had no effect. So, how can I fix the pressure and temperature for the boundary condition?

[dlparkhurst]

I'm going to say that is your problem, not PhreeqcRM. PhreeqcRM simply calculates reactions for the cells of the model at the specified temperature and pressures. If you set the temperature and pressure for the last cell of the model, those are the values that PhreeqcRM will use. If you plot temperature or pressure, they will simply be the values that you specified (excepting perhaps fixed volume GAS_PHASE calculations).

Your transport code implements the boundary conditions for the flow and transport equations.

[GeeqC]

Please let me rephrase the question:

I am not sure if the pressure and temperature is set correctly in Phreeqc for the boundary conditions.
Solution 2 is used for the boundary conditions as follows:

std::vector<double> bc_conc;
std::vector<int> bc2;
std::vector<double> bc_species;
bc2.push_back(2);
status = phreeqc_rm.InitialPhreeqc2SpeciesConcentrations(bc_conc, bc2);

The question is, whether the pressure and temperature are updated for Solution 2. From the results it seems that this is not the case. Whatever value I write for pressure and temperature, the outcome is always the same. Is there anything else I must do, so that Solution 2 with the modified temperature and pressure is used?

[dlparkhurst]

OK, so solution 2 is defined in the initial_phreeqc instance, and you are transferring it to PhreeqcRM. However, once it gets to PhreeqcRM, unless you put it in one of the cells, it is not touched by PhreeqcRM; the concentrations are simply for your boundary conditions.

However, solution 2 does not change just because you specify REACTION_TEMPERATURE or REACTION_PRESSURE.

So, if you you want the the same element concentrations, but only the temperature and pressure changed, you can use SOLUTION_MODIFY.

If you want a solution in equilibrium with some equilibrium phases, you will have to run a calculation:

Code: [Select]

USE solution 2
USE equilibrium_phases 2
REACTION_TEMPERATURE
100
REACTION_PRESSURE
50
SAVE solution 2
END


[GeeqC]

Thanks! This was very helpful.
The problem was: I did not SAVE solution 2.

Now my comment to the mineral reactions:

Example 6 in the manual suggests three different ways to solve the problem. If I understand this correctly, the first two methods are thermodynamic, using EQUILIBRIUM_PHASES. The third method is kinetic (using KINETICS). If I want to simulate a system with minerals that are far from equilibrium, I need to use the third method. Here, however, the problem is that Al becomes limiting, due to its low concentration. Even if Al consumed by one mineral is balanced by dissolution of another mineral, the reaction would come to a halt, due to numerical reasons. The time step would have to be probably thousand times smaller, to simulate this process realistically.

As a workaround I have now coupled the mineral reactions, so that the Al production and consumption is balanced. The reactions are simple source and sink terms in the transport equation.

[dlparkhurst]

I don't know if it is worth it, but to be a little more efficient, you could create a second PhreeqcRM for the boundary conditions. You could set temperature and pressure, RunCells, and then transfer the boundary conditions for the transport calculations.

Kinetics with small concentrations of Al may or may not be a problem. CVODE is supposed to handle stiff sets of rate equations. You won't know unless you try it.

[GeeqC]

Thanks for the tip. I will consider this for a future project (as it would mean that I have to change entirely the transport code).

As an additional test of the present model, I would like to re-equlibrate the final array of solutions with the atmosphere at 1 atm pressure. This would mimick the sampling procedure where probably CO2 is outgassing, which changes the DIC concentration and pH. Is there a way that the temperature and pressure can be changed and a gas phase can be added at the end of the run (even if the equilibrium and gas phases in the initial solution are turned off)?

ic1 = 1;            // Solution 1
ic1[N + i] = -1;       // Equilibrium phases none

[dlparkhurst]

I don't know how your transport model is set up, but if you are getting boundary conditions through the PhreeqcRM intial phreeqc instance, I thought it might be easier to change temperature and pressure if you put them in a second PhreeqcRM and then you could set temperature and pressure easily.

If you have n cells in PhreeqcRM, you can run the string "EQUILIBRIUM_PHASES 1-n; CO2(g) -3.4 10; END" for all workers, and the next RUN_CELLS will equilibrate with atmospheric CO2. For a parallel run, you will be defining an EQUILIBRIUM_PHASES for each cell in each thread or process, but it is the last step and everything will disappear at the end of the run. You could get fancier and add EQUILIBRIUM_PHASES appropriate for the cells in each thread or process, but for this case, I don't think it is worth the effort.

You can delete all the gas phases in all the cells if you RunString with the string "DELETE; -gas_phase; END

Temperature and pressure are set with SetTemperature and SetPressure.

[GeeqC]

Yes, we have n cells in PhreeqcRM and would like to run the string "EQUILIBRIUM_PHASES 1-n; CO2(g) -3.4 10; END" for all workers. How can this be done? Is it necessary to define a new PhreeqcRM instance and a new initial solution? And how does the program know to apply the equilibration to the array of the previously calculated solutions?

Is there perhaps a similar case in the example files? We would be very thankful for further guidance.

[dlparkhurst]

Just set a string equal to "EQUILIBRIUM_PHASES 1-n; CO2(g) -3.4 10; END" and use RunString for the workers.

Code: [Select]

std::string str = "EQUILIBRIUM_PHASES 1-n; CO2(g) -3.4 10; END"
phreeqc_rm.RunString(true,false,false,str);
phreeqc_rm.RunCells();

The definition of EQUILIBRIUM_PHASES will define EQUILIBRIUM_PHASES 1-n redundantly in every thread or process, so if you have 8 threads in OpenMP, you will have 8 copies of each EQUILIBRIUM_PHASE user number. However, RunCells will only run the cells for which a thread has SOLUTION definitions (roughly 1/8 of the cells in this example), and since this is the last calculation, you don't have to be too careful with efficiency.

[GeeqC]

This runs fine now. Thanks!

Now another error was found. We set the temperature and pressure gradients using:
        phreeqc_rm.SetTemperature(temperature);
        phreeqc_rm.SetPressure(pressure);

However, in the second cell from the top in our array, the P and T of the lower boundary condition were set, instead of the values from the gradients.

May this problem relate to your recommendation above, to use a new instance for P and T? But this part I have not quite understood.

[dlparkhurst]

If you are not using a mapping between your model cells and the chemistry cells, then the Workers have the same number of cells as the transport model, although split among multiple Workers if you are running in parallel. In the 1-1 mapping case, SetTemperature will assign the each temperature of the input array to the corresponding chemistry cell. SetPressure is the same, plus the pressure of the gas phase is also set for fixed_pressure gases.

If you do have a mapping, then there may be more cells in the transport model than there are chemistry cells, but the the temperatures and pressures are still set according to the mapping.

I don't know what has caused your problem. My guess is that you are mixing up the InitalPhreeqc instance and Worker instances. Did you RunString on the Workers that reset the temperature and pressure of cell 2? I don't see how SetTemperature and SetPressure could produce anything other than the values of the input arrays.

As for boundary conditions, you probably don't want to use another PhreeqcRM if you haven't gotten the first one working correctly. The advantage would be simplicity of determining determining the composition of a boundary solution that needs to be in equilibrium with other reactants at varying T & P.

[GeeqC]

In the meantime we could localize the problem:

The values of the lower boundary condition appear in the second cell from the top. This seems related to the way how we set the pressure and temperature for the boundary condition (solution 2):

        oss << "USE solution 2 " << "\n";
        oss << "REACTION_PRESSURE" << "\n";
        oss << pressure[N-1] << "\n";
        oss << "REACTION_TEMPERATURE" << "\n";
        oss << temperature[N-1] << "\n";
        oss << "SAVE solution 2 " << "\n";
        oss << "END\n";

If these lines are removed, the correct values appear in the second cell. But then the boundary condition solution is not speciated at the correct pressure and temperature.

[dlparkhurst]

If you are not running parallel, which I assume you are not for development, then if you run the script that you show through the Worker (RunString(true, ...)), then yes, model cell 2 will be set to the temperature and pressure that you specify in the script.

For all your model cells, SetTemperature and SetPressure can be used to set the temperature and pressure of the calculation.

So, if SetTemperature and SetPressure are not the answer, you are trying to set the temperature and pressure of a solution that is not one of your model cells? In that case, you probably want to do the calculation in either the Initial Phreeqc instance or the Utility instance. Typically, I would have the boundary condition solution to be used for the Transport calculation in the Initial Phreeqc instance, say solution 2 (this is not model cell 2, which is in the worker instances, it is just a temporary solution for calculations). You could then run your script (RunString(false, true, ...)) to determine solution 2 at the temperature and pressure that you want, and then get the concentrations for your transport model with InitialPhreeqc2Concentrations (or InitialPhreeqc2SpeciesConcentrations, if you are transporting species).

The element concentrations will not change with a change of temperature and pressure alone, but the species distribution will. If you are also equilibrating with a mineral or other reactant, then the element concentrations would also change with temperature and pressure. But however you do it, you can get the concentrations necessary for the boundary conditions of the transport model.

So, all of this would be to get concentrations to be used as boundary conditions in your transport model. After transport, you should have all the concentrations for the cells of the transport model, which are then sent to PhreeqcRM with SetConcentrations or SpeciesConcentrations2Module. And you can use SetTemperature and SetPressure to set the temperature and pressure of the RunCells calculation.

So, it comes down to the fact that I think you are not keeping the cells in the worker(s) separate from your calculation for the boundary condition. The calculation as you have described it is redefining the second cell of the model. Don't do this calculation in the worker instance(s).

[GeeqC]

The same problem as with solution 2 also occurs with solution 1. The initial solution 1 is set in the first (top) worker cell. This means that essentially the upper boundary conditions are set twice, one time in the reaction-transport model and one time in Phreeqc. Is there a way to avoid this duplication?

Also, in the equilibration with the atmosphere, at the end of the simulation, the first cell cannot be changed with:

oss << "EQUILIBRIUM_PHASES 1-" << N << "; CO2(g) -3.4 10";
oss << "END" << "\n";
phreeqc_rm.RunString(true,true,true,str);

How would it be possible to re-equilbrate the first cell?

[dlparkhurst]

Of course there is a way to avoid duplication, but I give up trying to explain it to you.

I don't think you understand the difference between the IPhreeqc instances: workers, initial phreeqc, and utility. When you do this

Code: [Select]

phreeqc_rm.RunString(true,true,true,str);

You are defining EQUILIBRIUM_PHASES 1-N for the workers--each worker has N equilibrium phases definitions, the initial iphreeqc instance--it has N equilibrium phases definitions, and the utility instance--it has N equilibrium phases definitions. So you have defined EQUILIBRIUM_PHASES 1 w+2 times, where w is the number of workers. When you use PhreeqcRM:RunCells, only one of these is used for each cell; it is the worker assigned a given cell number.

One last time. Do this:

(1) Use the initial phreeqc instance to calculate your boundary condition. Either define the conditions in a Phreeqc input file and RunFile(false, true, false, file), or some other way with RunString(false, true, false, str).

(2) use InitialPhreeqc2Concentrations or InitialPhreeqc2SpeciesConcentrations to get the boundary condition for your transport model.

(3) Run your transport step.

(4) Use SetConcentrations or SpeciesConcentrations2Module to set the concentrations in all of the cells in the worker instances (the cells of the model).

(5) If you want to equilibrate with CO2 in the final step, use your code, modified to so that only the workers run this string

Code: [Select]

oss << "EQUILIBRIUM_PHASES 1-" << N << "; CO2(g) -3.4 10";
oss << "END" << "\n";
phreeqc_rm.RunString(true,false,false,str);

(As explained above, it is inefficient to define a complete set of equilibrium phases for each worker because only one worker will use the definition for a cell. However, for the last calculation of the simulation, it should be adequate.)

(6) Use SetTemperature and SetPressure to set the temperature and pressure of the cells of the model (in the worker instances).

(7) Use PhreeqcRM:RunCells to do the calculation with the specified concentrations, equilibrium_phases, temperature, and pressure.

Make it so.

[GeeqC]

Sorry, I did not explain the problem very well. The simulation is a bit more complex, because there are initial conditions and boundary conditions. Here an overview:

     oss << "SOLUTION 1" << "\n"; //Initial conditions for the workers
     phreeqc_rm.RunString(true, false, false, str);

     oss << "SOLUTION 2" << "\n"; //Boundary condition for the transport model
     phreeqc_rm.RunString(false, true, false, str);

          R-T loop

     oss << "EQUILIBRIUM_PHASES 1-" << N << "; CO2(g) -3.4 10; END ";
     phreeqc_rm.RunString(true,false,false,str);

The problem is that the initial solution 1 is filled into the top worker cell and is not changed during equilibration with CO2.

[dlparkhurst]

The way it is designed to work is to define enough solutions (and other chemical reactants) in the Initial Phreeqc instance to be able to define all the initial conditions for the model , and then use InitialPhreeqc2Module to populate the model cells (in the workers). There are a few reasons for this plan. You can define just a few solutions in the initial phreeqc instance, you can define linearly varying conditions for the model cells with InitialPhreeqc2Module, but most of all, you don't have to know which worker contains which model cells when running parallel. InitialPhreeqc2Module assigns the distribution of model cells to workers.

You can use RunString(true, false, false) to assign initial conditions, but, when running parallel, you are either specifying solutions and reactants to workers that will not be used in calculations, or you have to do more work to figure out which cells are in which workers. Further, PhreeqcRM has load-balancing capabilities so that the distribution of cells among workers can change during the calculation.

Here is how I interpret what you have said.

Code: [Select]

     oss << "SOLUTION 1" << "\n"; //Initial conditions for the workers
     phreeqc_rm.RunString(true, false, false, str);

This will redefine the solution in the first cell of the model, whether you use SAVE or not. If you define only SOLUTION, then it will be the solution composition. If you define SOLUTION and also run reactions, then the reaction solution will be used in the model cell if you include SAVE solution 1, otherwise, the SOLUTION 1 definition will be used.

Code: [Select]

     oss << "SOLUTION 2" << "\n"; //Boundary condition for the transport model
     phreeqc_rm.RunString(false, true, false, str);

This will define a solution in the initial phreeqc instance, and you can use InitialPhreeqc2Concentrations to extract the concentrations for the reaction-transport model.

Code: [Select]

          R-T loop

OK, reaction transport is now done. Here, I expect that you would use SetConcentrations to specify the concentrations in the model cells before reactions. You could also set the temperature and pressure. SetConcentrations will overwrite SOLUTION 1, so, if SetConcentrations does not define the solution you want in cell 1, then you could define SOLUTION 1 after SetConcentrations instead of at the beginning as you have it. However, it would be reasonable to use the results of the R-T calculation.

Code: [Select]

     oss << "EQUILIBRIUM_PHASES 1-" << N << "; CO2(g) -3.4 10; END ";
     phreeqc_rm.RunString(true,false,false,str);

This last code will define equilibrium phases 1-N in every worker, but nothing will be calculated until you use RunCells. A better alternative is to define EQUILIBRIUM_PHASES 1 in the initial phreeqc instance, and then use InitialPhreeqc2Module to populate the cells. If you use "-1" for everything other than EQUILIBRIUM_PHASES, then nothing but EQUILIBRIUM_PHASES will change in the model cells. Using InitialPhreeqc2Module approach avoids duplicating definitions or having to know the distribution of cells among the workers.

You do not show RunCells, but, if you define EQUILIBRIUM_PHASES as you have done and then invoke RunCells, then whichever solution you have for cell 1 will be equilibrated with atmospheric CO2. So, I do not understand your statement

"The problem is that the initial solution 1 is filled into the top worker cell and is not changed during equilibration with CO2."

If cell 1 has some saturation, then it will be equilibrated with CO2(g) (and no other minerals) in the RunCells calculation.

[GeeqC]

Thank you for the detailed instructions. This was successful.

In the end, I actually used GAS_PHASE instead of EQUILIBRIUM_PHASES, so that the solutions can be equilibrated with a defined volume of headspace. The question is then, how much solution is equilbrated with the headspace? It makes a difference if I put 1 mL or 1 L of solution in contact with 1 L of headspace.

Also, I noticed that after equilibration all concentrations were a bit lower. Could this have something to do with the way we corrected the density of the solution (as discussed earlier, we used mmol/L units and the "-density calc" function)?

[dlparkhurst]

If you define GAS_PHASE in the initial IPhreeqc instance and populate the cells of the module with InitialPhreeqc2Module, then the gas phase is scaled when it is transferred to the reaction model with according to the setting of SetUnitsGasPhase in PhreeqcRM.

If you define GAS_PHASE directly in the workers, then you must consider the same factors that are discussed in the documentation for SetUnitsGasPhase.

[dlparkhurst]

As for the concentration change, you will have transfer of moles of elements between the gas phase and the solution. If you included H2O(g), then there may be an obvious change in the amount of water in the aqueous phase, but other gases could also affect the amount of water through hydrolysis reactions.

[GeeqC]

This discussion was very helpful to build our model linked to PhreeqRM, which is now published in:

Meister, P., Herda, G., Petrishcheva, E., Gier, S., Dickens, G.R., Bauer, C., and Liu, B. (2022) Microbial alkalinity production and silicate alteration in methane charged marine sediments: implications for porewater chemistry and diagenetic carbonate formation. Front. Earth Sci. 9, 756591.

https://doi.org/10.3389/feart.2021.756591