Example of using EQUILIBRIUM_PHASES_MODIFY in gas solubility problem

[Christinali91]

Hi David,

I am working on a reactive transport test case using phreeqc. My goal is to update CO2 concentration with EQUILIBRIUM_PHASES_MODIFY in every iteration. However, I am not sure how to correctly use EQUILIBRIUM_PHASES_MODIFY and didn't find a proper example online. Could you help me with it? Below is my code not working very well.


Thanks a lot for your help![/size]

Code: [Select]

SOLUTION 1  CaCl2
        -pressure 1000
        units            mmol/kgw
        temp             25.0
        pH               7.0     charge
        pe               12.5    O2(g)   -0.68
        Ca               0.6
        Cl               1.2
        C                0.001


EQUILIBRIUM_PHASES 1
  CO2(g) 0.1 # 1000 atm
END




EQUILIBRIUM_PHASES_MODIFY 1
 -component CO2(g)
 -si -2
 -moles 10
END


USE solution 1
USE EQUILIBRIUM_PHASES 1


SELECTED_OUTPUT
   -reset false
   -totals C
   -saturation_indices CO2(g)
   -file co2.dat
   USER_PUNCH
   -headings PRESSURE Fugacity PartialP CO2_g
   10 PUNCH  PRESSURE
   20 PUNCH  10^SI("CO2(g)")
   30 PUNCH  PR_P("CO2(g)")
   40 PUNCH  gas_vm
END


USE solution 1
USE EQUILIBRIUM_PHASES 1
EQUILIBRIUM_PHASES_MODIFY 1


SELECTED_OUTPUT
   -reset false
   -totals C
   -saturation_indices CO2(g)
   -file co2_2.dat
   USER_PUNCH
   -headings PRESSURE Fugacity PartialP CO2_g
   10 PUNCH  PRESSURE
   20 PUNCH  10^SI("CO2(g)")
   30 PUNCH  PR_P("CO2(g)")
   40 PUNCH  gas_vm
END







[dlparkhurst]

The EQUILIBRIUM_PHASES_MODIFY data block is described in the appendix of the manual. If you use DUMP to print out the raw form of an equilibrium phases definition, it is possible to change any value with EQUILIBRIUM_PHASES_MODIFY.

It was intended that only the values in the "candidate identifiers" section would be meaningful, and the rest would be work space variables that would be overwritten in the process of making the calculations. However, it appears that a "feature" was introduced when implementing Peng-Robinson for gases. You should specify both -si and -si_org to ensure that a new target saturation index (log partial pressure for a gas) is installed. -si is in the candidate list, but -si_org is not. For an ideal gas (wateq4f.dat for example only implements ideal gases), you need only specify -si, but for phreeqc.dat, Amm.dat, or pitzer.dat, it is necessary to specify both -si and -si_org if you want to change the target partial pressure of a gas.

The most likely item to change is the moles of the phase or the target si. The following changes both (for phreeqc.dat) as an example. By default, dump output goes to the file dump.out.

Code: [Select]

SOLUTION 1  CaCl2
        -pressure 1000
        units            mmol/kgw
        temp             25.0
        pH               7.0     charge
        pe               12.5    O2(g)   -0.68
        Ca               0.6
        Cl               1.2
        C                0.001
END
USER_PRINT
10 PRINT "P(CO2):        ", 10^SI("CO2(g)") / PR_PHI("CO2(g)")
20 PRINT "Fugacity(CO2): ", 10^SI("CO2(g)")
30 PRINT "Phi(CO2):      ", PR_PHI("CO2(g)")
USE solution 1
EQUILIBRIUM_PHASES 1
  CO2(g) 1.0 0.1          # 10 atm target pressure, 0.1 moles available to react
SAVE equilibrium_phases 1 # 0.01916 moles remain
DUMP
-equilibrium_phases 1
END
EQUILIBRIUM_PHASES_MODIFY 1
-component CO2(g)
   -moles   10   # moles increased from 0.01916 to 1.0
   -si      2
   -si_org  2   # target partial pressure 100 atm
END
DUMP
-append
-equilibrium_phases 1
END
USE solution 1
USE equilibrium_phases 1
SAVE equilibrium_phases 1
END
DUMP
-append
-equilibrium_phases 1
END


[Christinali91]

Hi David,

Thanks so much for your quick response. Now I get my phreeqc working. I have a follow-up question about EQUILIBRIUM_PHASES_MODIFY. I try to couple phreeqc with the PhreeqcRM interface. However, my EQUILIBRIUM_PHASES_MODIFY seems not working in the code.

In the initialization step, I have Ierr = RM_RunFile(Id_phreeqc,  1, 1, 1, 'IPhreeqc_in.txt') to read Phreeqc input as below:

Code: [Select]

SOLUTION 1-100
-water 1.0 ## water mass in each grid block
units mmol/kgw
temp 31.05
pH 7.0 charge
pe 12.5 O2(g) -0.68
Ca               0.6
Cl               1.2
C                0.001

EQUILIBRIUM_PHASES 1-100
  CO2(g) 1.0 0.1
END
Then I use EQUILIBRIUM_PHASES_MODIFY to modify -si and -si_org as below:

Code: [Select]

            WRITE(DUMMY_STR,*) K
            input_IHC_temp = 'EQUILIBRIUM_PHASES_MODIFY  '//TRIM(ADJUSTL(DUMMY_STR))//new_line(cc)
            !
            DO II=1, NO_HC_GC
                WRITE(DUMMY_STR,*)    2.0
                WRITE(DUMMY_STR_2,*)  10.0
                !WRITE(*,*) 'si, moles',DUMMY_STR, DUMMY_STR_2
                input_IHC_temp = input_IHC_temp //'-component  '//TRIM(ADJUSTL(NAME_HC_GC(II)))//new_line(cc)
                input_IHC_temp = input_IHC_temp // '-si  '//TRIM(ADJUSTL(DUMMY_STR))//new_line(cc)
                input_IHC_temp = input_IHC_temp // '-moles  ' //TRIM(ADJUSTL(DUMMY_STR_2))//new_line(cc)
                input_IHC_temp = input_IHC_temp // '-si_org '//TRIM(ADJUSTL(DUMMY_STR))//new_line(cc)
                !input_IHC = input_IHC //'END' //new_line(cc)
            ENDDO  !! for NO_HC_GC
            input_IHC = input_IHC // input_IHC_temp // new_line(cc)
        ENDDO
        input_IHC = input_IHC //'END' //new_line(cc)
        Ierr = RM_RunString(Id_phreeqc, 1, 1, 0, input_IHC)

I also print out the contents of this string as below. The output seems correct to me.

Code: [Select]

EQUILIBRIUM_PHASES_MODIFY  1
-component  CO2(g)
-si  2.000000
-moles  10.00000
-si_org 2.000000

EQUILIBRIUM_PHASES_MODIFY  2
-component  CO2(g)
-si  2.000000
-moles  10.00000
-si_org 2.000000
...
However, when I use SELECTED_OUTPUT and USER_PUNCH to get CO2_SI, the value changes but doesn't meet my target SI.

      Selected output:
            1                     C:      1.1020E-01
            2                    Ca:      5.9825E-04
            3                    Cl:      1.1965E-03
            4             TOT("Ca"):      6.0000E-04
            5             TOT("Cl"):      1.2000E-03
            6           TOT("C(4)"):      1.1052E-01
            7          DELTA_CO2(g):      0.0000E+00
            8             SI_CO2(g):      5.7652E-01
            9                Charge:     -7.6349E-18
           10              Excess_O:      4.7883E-04
           11              Excess_H:      2.2112E-01
           12                   H2O:      5.5344E+01

Could you give me some hints about it?

Thank you so much for the time and help!


Best,
Christina

[dlparkhurst]

It's tough to know what goes on between transport and chemistry. All I can say is that this PHREEQC file shows what I think should happen if you do the equivalent in PhreeqcRM using the phreeqc.dat database. Note that RUN_CELLS does a SAVE for both solution and equilibrium phases.

Note that 0.1 mol of CO2 is not sufficient to generate a partial pressure of 10 in the first RUN_CELLS. In the second RUN_CELLS, 10 moles is sufficient to produce a partial pressure of 100 atm.

Note there is an inconsistency in PHREEQC that evolved with the addition of Peng-Robinson. When you set -si in EQUILIBRIUM_PHASES and EQUILIBRIUM_PHASES_MODIFY, you are setting the log10 partial pressure. The Basic function SI returns the log10 fugacity. Fugacity and pressure are equal when using the ideal gas law, which was implemented in most PHREEQC databases. However, fugacity and pressure are not equal with Peng-Robinson databases (phreeqc.dat, Amm.dat, and pitzer.dat). With these databases, the Basic function SI("CO2(g)") will not equal the value for -si in EQUILIBRIUM_PHASES_MODIFY. To get the value of -si in EQUILIBRIUM_PHASES_MODIFY, you need to calculate it from fugacity and phi, for example, logP = SI("CO2(g)") - log10(PR_PHI("CO2(g)")).

Code: [Select]

USER_PRINT
10 PRINT "P(CO2): ", 10^SI("CO2(g)") / PR_PHI("CO2(g)")
SOLUTION 1
-water 1.0 ## water mass in each grid block
units mmol/kgw
temp 31.05
pH 7.0 charge
pe 12.5 O2(g) -0.68
Ca               0.6
Cl               1.2
C                0.001
END
EQUILIBRIUM_PHASES 1
  CO2(g) 1.0 0.1
END
RUN_CELLS
-cell 1
END
EQUILIBRIUM_PHASES_MODIFY  1
-component  CO2(g)
-si  2.000000
-moles  10.00000
-si_org 2.000000
END
RUN_CELLS
-cell 1
END


[Christinali91]

Hi David,

Thanks so much for your explaination. I did miss the difference between Fugacity and pressure in Peng-Robinson. Now the fugacity value seems correct to me. I also want to get the mole or mole fraction of CO2 using SELECTED_OUTPUT and USER_PUNCH.

I use TOT("C(4)") for CO2 in water. As for CO2 in gas phase, I tried GAS(CO2(g)) and EQUI("CO2"). However, GAS(CO2(g)) seems only working for GAS_PHASE but not EQUILIBRIUM_PHASES while the value of EQUI("CO2") is always zero in my code. Is there other keywords I can use to get mole of CO2(g) during the simulation?


Thanks again!

Sincerely,
Christina

[dlparkhurst]

Using TOT("C(4)") will give you the concentration of all carbonate species (including CO2) in units of mol/kgw. The total number of moles would be TOT("C(4)") * TOT("water"), where TOT("water") is the mass of water (alternatively, TOTMOL("C (4)") would give the same quantity).

If you want the total moles of CO2(aq), you would use MOL("CO2") * TOT("water").

EQUI("CO2(g)") returns the number of moles of the phase CO2(g) in the current state of EQUILIBRIUM_PHASES.

[Christinali91]

Hi David,

Thanks for the timely reply. Now I can print out all the varables as below. However, I find that EQUI("CO2(g)") + aq_CO2(MOL("CO2") * TOT("water")) is larger than TOTMOL("C(4)"). Is this possible?  I thought TOTMOL("C(4)") should always be equal or larger than EQUI("CO2(g)") + aq_CO2. This is no transport involved at this point.

If EQUI("CO2(g)") is not the correct keyword to get moles of CO2 in the gas phase(CO2(g))? Is there any other keyword I can use?

Thanks a lot for your help!

Code: [Select]

Cell number            1
      Components:
           1         H2O:    54.40945791
           2           H:     0.62618144
           3           O:     0.54892148
           4      Charge:     0.00000000
           5           C:     0.43063331
           6          Ca:     0.00058927
           7          Cl:     0.00117855
      Selected output:
            1                     C:      4.3776E-01
            2                    Ca:      5.9902E-04
            3                    Cl:      1.1980E-03
            4       EQUI("CO2(aq)"):      0.0000E+00
            5        EQUI("CO2(g)"):      7.1460E-02
            6        TOTMOL("C(4)"):      2.7881E-01
            7       Fugacity_CO2(g):      9.4863E-01
            8             SI_CO2(g):      9.7710E-01
            9                 aq_CO2:      2.7505E-01
           10                Charge:      7.8943E-19
           11              Excess_O:      6.3654E-01
           12              Excess_H:      5.5800E-01
           13                   H2O:      5.5309E+01


[dlparkhurst]

TOT, TOTMOL, and MOL refer only to the aqueous solution, so EQUI values are not included in TOTMOL.

You can use SYS to find the total number of moles in all phases of a cell for an element [SYS("C")] or redox state [SYS("C(4)")].