How is the solution volume calculated in phreeqc?

[Wonwoo Yoon]

Hello.

I’m currently working on coupling COMSOL and phreeqc. Up until this point, the transition between molarity and molality posed no issues due to the system’s low concentration and isothermal assumption. Consequently, I simply assumed mol/kgw equaled mol/L. However, my current challenge involves dealing with seawater in a non-isothermal system, requiring a precise unit conversion between molarity and molality.

Presently, I convert mol/L in COMSOL to mol by multiplying the porosity and saturation. Then I incorporated these values into phreeqc through SOLUTION_MODIFY. After the calculation of phreeqc, the total mol in the system is obtained through SELECTED_OUTPUT (e.g., TOTMOL(“Ca”)). I tried to get the molarity by dividing the obtained total mol by the solution volume. However, discrepancies arise between COMSOL and phreeqc in water density calculations, which might result in different solution volumes. To maintain consistency, I aim to unify the density calculation.

To achieve this, I seek insight into phreeqc’s methodology for computing water density and solution volume. My assumption is that water density is determined using a specific equation of state (I don’t know which of them was used, though), while the solution volume results from the sum of water volume (by multiplying water mass and water density) and the volume effect of species. I would greatly appreciate specific information regarding this. Any additional discussions regarding the mentioned coupling sequence are also welcomed.

Thanks.

Wonwoo Yoon

[dlparkhurst]

You can find discussions of the density calculation in the SOLUTION_SPECIES description in the version 3 manual and in the Geochimica et Cosmochimica Acta paper "Equations for calculating hydrogeochemical reactions of minerals and gases such as CO2 at high pressures and temperatures", by Appelo, Parkhurst, and Postma.

You do not mention PhreeqcRM, which is the package that you should be using for reactive-transport calculations. It is designed to be the reaction module for multicomponent transport models. PhreeqcRM has been used successfully in many reactive transport applications, including several using COMSOL. It also has parallelization by either OpenMP or MPI.

PhreeqcRM does the conversions between the transport model units, typically mg/L or mass fraction, and Phreeqc units. The method GetConcentrations provides an array of concentrations in mg/L (for example), and SetConcentrations takes the transported concentrations in mg/L and updates the Phreeqc solutions. More details about the conversion are available at https://water.usgs.gov/water-resources/software/PHREEQC/documentation/phreeqcrm-html/index.html.

The Phreeqc-calculated densities are available with the PhreeqcRM method GetDensity. You can either lag the values or iterate between chemistry and transport in your numerical method to account for density changes in the transport equations.

[Wonwoo Yoon]

Thank you for your kind reply.

I will figure out that phreeqcRM can be utilized in our current system.

Several works including code validation by comparing results to that of other codes and some main results production have been already done. Therefore I have to decide whether to change the whole system to phreeqcRM, or just figure out how the phreeqcRM deals with such coupling problem, and apply it to our code.

Regardless of that, your reply was greatly helpful to me. Thanks again for that.

Sincerely, Wonwoo Yoon

[dlparkhurst]

Maybe not the answer you want to hear, but I think it would be worthwhile to switch to PhreeqcRM. PhreeqcRM has been used many times without significant bug reports, so the code should be reliable for what it does. It has proved sufficient for many types of reactive transport modeling, including density and multiphase applications. The coupling for calculations is fairly straightforward; most of your issues will be input and output. In the long run, it should be easier to justify your methods and maintain your code.

I will add that a new BMIPhreeqcRM version is in the works. It implements the Basic Model Interface, which adds functionality for coupling with other BMI codes, but also simplifies retrieving model results. All current capabilities of PhreeqcRM are unchanged.

[Wonwoo Yoon]

Great to hear that!

Actually, I installed the phreeqcRM and figured out how to run the C library in MATLAB (which is our code's main body), yesterday.

However, there are still many things to do, mainly about changing the method of importing and exporting the input and output from COMSOL (I think the current way of importing and exporting data of our code is not efficient when it comes to phreeqcRM.)

Again, thanks for your help!

Wonwoo Yoon