Use RM_SetConcentrationsand RM_GetConcentrations in two-phase problem

[Christinali91]

Hi David,

I am working on a problem simulating CO2 dissolved in water. I use PhreeqcRM for chemical reaction calculation to get species distribution including the amount of CO2 and H2O species left after chemical reaction and then use my own EOS solver for phase equilibrium calculation. I'll iterate between them to make sure the results are consistent. However, I am confused about the results I get from Phreeqc.

Since this is a two-phase problem, I use RM_SetConcentrations and RM_GetConcentrations to update all the dissolved C related species. To be more specific, the concent array I use in RM_SetConcentrations includes TOT_C, charge, TOT_H, TOT_O meanwhile I use EQUILIBRIUM_PHASES_MODIFY to update CO2 in the gas phase. Then I ran phreeqc to get results(as below).

Code: [Select]

GLOB_TIME           0
 Cell number            1
      Components:
           1           H:    98.07812508
           2           O:    56.86409841
           3      Charge:    -0.00000000
           4           C:     3.91251793
      Selected output:
            1                     C:      4.4155E+00
            2            MOL("CO2"):      3.7878E+00
            3         MOL("(CO2)2"):      3.1951E-01
            4            MOL("H2O"):      5.5506E+01
            5        EQUI("CO2(g)"):      0.0000E+00
            6       10^SI("CO2(g)"):      2.9696E+02
            7                Charge:     -2.9063E-14
            8               Total_H:      1.1069E+02
            9               Total_O:      6.4175E+01
There are a few things I don't quite understand:
(1) The values of H, O and C from RM_GetConcentrations are smaller than the results I got from SELECTED_OUTPUT. In SELECTED_OUTPUT, I use the PUNCH keyword as below. I expected the results were consistent with the results from RM_GetConcentrations, but they were not.

Code: [Select]

          10                      PUNCH TOTMOLE("C")
          80                      PUNCH TOTMOLE("H")
          90                      PUNCH TOTMOLE("O")
(2) The value of dissolved CO2(MOL("CO2")+2MOL("(CO2)2")) is larger than TOTMOLE("C"), which I don't think is correct. TOTMOLE("C") should be larger than the value of dissolved CO2 since there are other carbonates and bicarbonate species. I also tried [SYS("C(4)")], the value of [SYS("C(4)")] is equal to TOTMOLE("C").

Code: [Select]

USER_PUNCH  1
-Heading   C                                                           
-Heading   MOL("CO2")                                                 
-Heading   MOL("(CO2)2")                                               
-Heading   MOL("H2O")                                                 
-Heading   EQUI("CO2(g)")                                             
-Heading   10^SI("CO2(g)")                                             
-Heading   Charge
-Heading   Total_H
-Heading   Total_O
          10                      PUNCH TOTMOLE("C")
          20                      PUNCH MOL("CO2")
          30                      PUNCH MOL("(CO2)2")
          40                      PUNCH MOL("H2O")
          50                      PUNCH EQUI("CO2(g)")
          60                      PUNCH 10^SI("CO2(g)")
          70                      PUNCH CHARGE_BALANCE
          80                      PUNCH TOTMOLE("H")
          90                      PUNCH TOTMOLE("O")
Could you give me some hints on these or maybe some suggestions on how I can check if my results are reasonable?

Thanks so much for your help!

Best,

Christina

[dlparkhurst]

There are several different units involved in the values you cite.

TOTMOL returns the number of moles in solution.
SYS returns the number of moles in the cell including all reactants (EQUILIBRIUM_PHASES,   EXCHANGE, etc.)
TOT returns mol/kgw in solution.
MOL returns mol/kgw in solution.
GetConcentrations depends on the value given in SetUnitsSolution. Options are 1, mg/L; 2 mol/L; or 3, mass fraction, kg/kgs.

Some relavent Basic functions are

TOT("water") is the mass of water in kilograms
SOLN_VOL is the volume of the solution in liters (accurate for phreeqc.dat, Amm.dat, and pitzer.dat)
RHO is the density of the solution, kg/L.

TOTMOL = TOT * TOT("water")
For an aqueous species mol/L = MOL*TOT("water")/SOLN_VOL
For a total element, mol/L = TOTMOL/SOLN_VOL = TOT*TOT("water")/SOLN_VOL
Mass of solution, kg = SOLN_VOL*RHO

Conversions for GetConcentrations and SetConcentrations are described in documentation for SetUnitsSolution at https://water.usgs.gov/water-resources/software/PHREEQC/documentation/phreeqcrm/. Note that PhreeqcRM internally saves solutions as moles of elements and charge.

[Christinali91]

Hi David,

Thanks so much for your reply! I have a follow-up question about SOLN_VOL.

How does Phreeqc get SOLN_VOL? And how can I set the value of SOLN_VOL for my problem?

Thanks again!

Best,

[dlparkhurst]

If you really want to know, density is calculated as follows:

rho = rho_0 * (1e3 + M_T / mass_water_aq_x) / (rho_0 * V_solutes / mass_water_aq_x + 1e3);

where rho_0 is the density of pure water at the temperature and pressure, M_T is the mass of solutes in grams, mass_water_aq_x is the mass of water, V_solutes is the sum over all species of moles*Vm, and Vm is calculated from parameters in SOLUTION_SPECIES input.

For phreeqc.dat, Amm.dat, and pitzer.dat, Vm is defined for each species. For other databases, Vm is not defined (0.0), so the V_solutes=0.

Solution volume is then calculated as:

vol = 1e-3 * total_mass / rho;

where total_mass is again the total mass of solutes in grams.

PhreeqcRM gives options for calculating the conversions in SetConcentrations and GetConcentrations. See SetUnitsSolution for details. If you want to avoid the issues with solution volume and density, maybe you should use the mass fraction option.