Exporting Equilibrium Phases moles in Reactive Transport with PhreeqcRM

[Sara Tabrizi]

Hello all,
I have coupled PhreeqcRM with a transport code for a test case containing Calcite and Dolomite dissolution/precipitation in the equilibrium phase. I would like to export the moles of equilibrium phases (Calcite and Dolomite) for model verification. Although it is very easy to export the moles per volume for equilibrium phases in the batch Phreeqc file by using the command "-equilibrium_phases" in the selected output data block, I cannot find the command to export these values in PhreeqcRM. Even though there are classes such as "GetEquilibriumPhasesCount" or "GetEquilibriumPhases" to get the numbers and names of the phases, I cannot get the values of moles for each equilibrium phase. I have been studying different classes in PhreeqcRM to find the relevant command for a while without success. I appreciate it if anyone who knows lets me know about exporting the moles of equilibrium phases during the reactive transport.
Thank you in advance,
Sara Tabrizinejadas

[dlparkhurst]

In a few months, it should be much easier. We are working on a BMI version of PhreeqcRM that will provide a lot more variables through the BMI function GetValue.

For now, you will have to do it the old fashioned way by defining a selected output file. You will need to define SELECTED_OUTPUT n and USER_PUNCH n (if needed) data blocks and run them in the worker instances. There is the PhreeqcRM method SetSelectedOutputOn true/false that you can switch on for time steps you want to get the selected-output data.

Finally, there are methods to extract the selected-output data: GetSelectedOutputColumnCount, GetSelectedOuputHeadings, GetSelectedOutputRowCount, and GetSelectedOutput that allow you to retrieve and process the selected-output data.

You noted the lists of names of EQUILIBRIUM_PHASES, exchange species, surface species, and others are available, so you could automatically generate output by processing these names and writing a SELECTED_OUTPUT definition that you then process with the workers. The file advection_cpp.cpp in the Tests directory of the PhreeqcRM distribution has a rough version of how to generate a SELECTED_OUTPUT definition automatically in the subroutine example_selected_output.

[Sara Tabrizi]

Dear Dr.Parkhurst,
Thank you very much for your reply. I have been trying to obtain the values of equilibrium phase moles with this approach. To do so I have defined equilibrium_phases in the Selected_output data block of the input .pqi file and this file has been run in the worker instances. I have used the following lines at the end of .pqi file:
"
SELECTED_OUTPUT 1-500
-distance true
-time true
-equilibrium_phases Calcite Dolomite
"
Also in the transport code I use the following lines to define the selected output:
"   
IRM_RESULT status;
status = phreeqc_rm->SetSelectedOutputOn(1);
"

But apparently, this is not enough to define the selected output in the Phreeqcrm part of the code. Because when I call the class "GetSelectedOutputColumnCount()", I receive the error that the Selected output is not found. Could you please let me know what has been done wrong or what are the other necessary definitions to obtain the desired selected output?
Thank you very much,
Sara

[dlparkhurst]

My bad. I forgot to explain that you need to select the selected-out definition that you wan

You can define multiple SELECTED_OUTPUT data blocks, which can be processed with PhreeqcRM. First, you should use only one user number for each (1, not 1-500). The SELECTED_OUTPUT/USER_PUNCH definitions must be run by the workers.

Here is the missing step--with the PhreeqcRM method SetCurrentSelectedOutputUserNumber you choose the SELECTED_OUTPUT that you want to retrieve.

Here is a snippet of code from the html documentation (https://water.usgs.gov/water-resources/software/PHREEQC/documentation/phreeqcrm-html/index.html or in the distribution directory/Doxygen/html/index.html).

Code: [Select]

for (int isel = 0; isel < phreeqc_rm.GetSelectedOutputCount(); isel++)
{
  int n_user = phreeqc_rm.GetNthSelectedOutputUserNumber(isel);
  status = phreeqc_rm.SetCurrentSelectedOutputUserNumber(n_user);
  std::vector so;
  int col = phreeqc_rm.GetSelectedOutputColumnCount();
  status = phreeqc_rm.GetSelectedOutput(so);
  // Process results here
}


[Sara Tabrizi]

Thanks again for your reply. I think that this is not the problem. I changed what I had written in the .pqi file to the following lines :
"
SELECTED_OUTPUT 1
-file RUN1.sel
-distance true
-time true
-equilibrium_phases Calcite Dolomite
"
Also, I call the following required classes before running the .pqi file in the worker instance :
SetSelectedOutputOn(true);
SetCurrentSelectedOutputUserNumber(1);

However, it appears to me that the output of  "GetSelectedOutputCount()" is zero and that is the reason why the selected output is not defined in the code. I do not understand what I am doing wrong, I appreciate it if you could give me a hint if something wrong in my procedure is pointing out to you.
Sorry for the repetitive question and thanks for your patience,
Sara

[dlparkhurst]

A couple of things first, the -file definition does not apply during time-stepping with PhreeqcRM. You will have to write selected-output results to a file. Also, the distance applies in the initial Phreeqc run (if you were using TRANSPORT, not likely), but will not be meaningful during time stepping. I'm not sure if -time will be correct or not.

Have you run the pqi file so that the SELECTED_OUTPUT definition is defined to the workers? I don't think that you have, or else the count would not be zero.

RunFile(true, true, false, filename)

[Sara Tabrizi]

I have changed the SELECTED_OUTPUT datablock by removing the -file, -time and -distance. Also I believe that I am running the .pqi file within the command RunFile(true, true, false, filename). I am using the following order:
"
status = phreeqc_rm->SetSelectedOutputOn(true);
status = phreeqc_rm->RunFile(true, true, false, "filename");
status = phreeqc_rm->SetSelectedOutputOn(true);
status = phreeqc_rm->SetCurrentSelectedOutputUserNumber(1);
int a= phreeqc_rm->GetSelectedOutputCount();
"
However, the output of the last line (int a) is zero. I think it might be due to writing the selected output results to a file. without the command -file in the selected_output data block what would be the way to do this?

[dlparkhurst]

Have you run RunCells before you are trying to extract the data?

Follow the sequence of method calls in the Advection example in Tests/Advection_cpp.cpp or Tests/Advection_f90.f90.

[Sara Tabrizi]

The output of the selected output count is not zero after Runcells. I should have called them after Runcells. Sorry about this mistake and thank you so much for your help.

[dlparkhurst]

There are dependencies in the order of method calls that probably need more documentation or warning messages. But, I'm not sure how to do it, and it probably won't get done.

[Sara Tabrizi]

I think that the documentation is very well written in your publication:
https://doi.org/10.1016/j.advwatres.2015.06.001

[Sara Tabrizi]

Hello again dear Dr.Parkhurst,
I exported the moles of equilibrium reaction with the classical method of Selected output data block. However, the exported values from Phreeqcrm are less than the values that are supposed to be. I think this is due to the several chemical reactions considered in Phreeqcrm that in the transport code, I have not considered their resulting species. That is why I have to modify the species concentrations and assign the values that are not considered in the transport code to the ones that are considered. That is working fine for the dissolved species. However, for the minerals, I still have different results than the batch Phreeqc result because they are exported based on a different chemical system. I was wondering what formulation has been used in Phreeqc for calculating the moles of equilibrium phases? that would be a great help if I know the formulation, so I would be able to calculate the moles of equilibrium phases based on the chemical system corresponding to the transport system.
Thank you in advance,
Sara

[dlparkhurst]

When using the same versions of Phreeqc and PhreeqcRM, the calculation engine should be the same.

You should use the same database when comparing Phreeqc and PhreeqcRM.

If you simply equilibrate the reactants with the solutions without transport (just RunCells before transport calculations), you should get the same results with Phreeqc batch calculation with the same solution and reaction definitions.

When you add flow and transport, you are getting into more possibilities for differences.

[dlparkhurst]

I'm not sure what you did, but I don't think it is a good idea. I don't think the gain in speed would be worth it, and although your reaction network may be simpler, comparing to PHREEQC calculations appear to be problematic. You are adding extra coding to do some sort of transformation, and probably eliminating more complicated types of reactions, like solid solutions, surfaces, and exchange without even more coding.

A number of moles of mineral is assigned in the EQUILIBRIUM_PHASES definition. Phreeqc dissolves or precipitates the mineral to the point that the saturation index is zero. The number of moles of mineral is changed, and the solution composition increases or decreases the concentrations of the elements from the mineral reaction. In the case of a more complicated system with more reactants, the elements added or removed affect the composition of all the reactants (surfaces, exchangers, solid-solutions, gas phase, and others).

Note that PhreeqcRM uses FindComponents to return only elements that are found in the Initial Phreeqc instance, usually only the elements used in the solutions and reactants in a PHREEQC input file need to be transported. If you do not include Fe in any solution or reactant (EXCHANGE, EQUILIBRIUM_PHASES, etc), then it will not be in the list provided by FindComponents of elements that need to be transported. So, if
you are careful with your Phreeqc input file, you will get a minimum set of components to transport.

If you are trying to remove H, and O from being transported, I think your approach will fail. Charge balance could be removed if you are very careful to use only charge-balanced solutions and you do not use some surface reactions that introduce a charge imbalance. I don't think the savings from removing one transport equation is worth the risk of unexpected results if the conditions of charge balance are not met.

Another approach, if you want a simpler reaction network, simplify the database. In that case, you can compare your calculations to Phreeqc directly, and you won't have the extra coding.