Resources for setting up input file for marine CO2 modeling?

[mcringham]

Hello, I'm an oceanographer new to Phreeqc, and I'm struggling to wade through examples in the user manual. I want to model what I think should be a simple problem:

I have a tank of seawater open to the atmosphere. I know my seawater's typical starting chemistry (Na, Mg, etc concentrations are ~average seawater values), temperature, salinity, Dissolved oxygen, pH, DIC, and total alkalinity. I also know the atmospheric pCO2.

I add a known amount of base (NaOH) to this seawater, raising the pH and total alkalinity, and impacting the inorganic carbon speciation. I leave the tank alone, and over time, pH decreases, total alkalinity remains stable (no observed precipitation to affect this), and DIC increases as CO2 from the atmosphere intrudes into the tank.

In the end, I would like to know the estimated pH and carbon species concentrations as the seawater comes to equilibrium with the atmosphere on a timescale of weeks-months. I have experimental lab data for all of this that I would like to validate with this Phreeqc simulation. In all of this, I am interested only in the inorganic carbon chemistry and am neglecting other reactions.

Where would you send an absolute beginner for this type of problem? Are there examples available that I can modify to this purpose?

[dlparkhurst]

Here is something you can work from. In this script, I have used a seawater composition with SOLUTION. This produces an initial solution calculation with a distribution of species. This solution is then equilibrated with atmospheric CO2 with EQUILIBRIUM_PHASES in a reaction calculation, which adjusts the TDIC and pH. This reacted solution is saved as solution 2. Solution 2 is reacted with 0.01 mol of NaOH in a closed system (no CO2(g) in EQUILIBRIUM_PHASES), with the result saved as solution 3. Finally, solution 3 is equilibrated with atmospheric CO2.

Code: [Select]

SOLUTION 1  SEAWATER FROM NORDSTROM AND OTHERS (1979)
        units   ppm
        pH      8.22
        pe      8.451
        density 1.023
        temp    25.0
        redox   O(0)/O(-2)
        Ca              412.3
        Mg              1291.8
        Na              10768.0
        K               399.1
        Fe              0.002
        Mn              0.0002  pe
        Si              4.28
        Cl              19353.0
        Alkalinity      141.682 as HCO3
        S(6)            2712.0
        N(5)            0.29    gfw   62.0
        N(-3)           0.03    as    NH4
        O(0)            1.0     O2(g) -0.7
EQUILIBRIUM_PHASES 1
    CO2(g)    -3.4 10
SAVE solution 2
END
USE solution 2
REACTION 1
    NaOH       1
    0.01 moles
SAVE solution 3
END
USE solution 3
EQUILIBRIUM_PHASES 1
    CO2(g)    -3.4 10
END


[mcringham]

Thank you! This was really valuable advice. I've edited my code and have some follow up questions. This is my starting code, where I am using a large tank of seawater, already ~equilibrated with the atmosphere, at a known starting pH (7.8514), DIC (2271.1 umol/kg), TA (2282.4 umol/kg), Temp, and Salinity. Air pCO2 above the tank is ~420ppm CO2. The other chemical values are estimates from your seawater example.

Code: [Select]

TITLE: Experiment 1 test
SOLUTION 1 Tank2 Seawater
    temp      19.34
    pH        7.8514
    pe        4
    redox     O(-2)/O(0)
    units     ppm
    density   1.024
    #Alkalinity 2282.4 ueq/kgs    as HCO3
    C(4)      2271.1 umol/kgs
    Ca        412.3
    Cl        19353
    Fe        0.002
    K         399.1
    Mg        1291.8
    Mn        0.0002 pe
    N(-3)     0.03 as NH4
    N(5)      0.29 gfw 62
    Na        10768
    O(0)      1 O2(g)      -0.7
    S(6)      2712
    Si        4.28
    -water    6328 # kg
END
USE solution 1
REACTION 1
    NaOH       1
    2.3225 moles
SAVE solution 2
END
USE solution 2
EQUILIBRIUM_PHASES 1
    CO2(g)    -3.5 10
END
I react this tank with a known amount of NaOH, then allow the tank to equilibrate with the atmosphere again. I know from measurements that the pH immediately after reacting with NaOH is ~8.36 (seawater scale), and pH after several weeks is ~8.07 and decreasing. I have pH, DIC, TA, etc throughout this experiment, so I'm trying to tune this code to approximate measurements.

Some questions:
1) I know pH, TA, and DIC. I choose to identify pH and DIC to calculate the rest of the system here, bypassing TA. But my output shows the following:

Code: [Select]

-------------------------------------------
Beginning of initial solution calculations.
-------------------------------------------

Initial solution 1. Tank2 Seawater

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

C(4)              2.354e-03   1.489e+01
Ca                1.066e-02   6.746e+01
Cl                5.657e-01   3.580e+03
Fe                3.711e-08   2.349e-04
K                 1.058e-02   6.694e+01
Mg                5.507e-02   3.485e+02
Mn                3.773e-09   2.387e-05
N(-3)             1.724e-06   1.091e-02
N(5)              4.847e-06   3.067e-02
Na                4.854e-01   3.072e+03
O(0)              4.833e-04   3.058e+00  Equilibrium with O2(g)
S(6)              2.926e-02   1.851e+02
Si                7.382e-05   4.672e-01

----------------------------Description of solution----------------------------

                                       pH  =   7.851   
                                       pe  =   4.000   
      Specific Conductance (µS/cm,  19°C)  = 46693
                          Density (g/cm³)  =   1.02482
                               Volume (L)  = 6399.25495
                        Activity of water  =   0.981
                 Ionic strength (mol/kgw)  =   6.769e-01
                       Mass of water (kg)  =   6.328e+03
                 Total alkalinity (eq/kg)  =   2.408e-03
                       Total CO2 (mol/kg)  =   2.354e-03
                         Temperature (°C)  =  19.34
                  Electrical balance (eq)  =   5.009e+00
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =   0.07
                               Iterations  =   7
                                  Total H  = 7.025028e+05
                                  Total O  = 3.520334e+05
------------------
End of simulation.
------------------

-----------------------------------------
Beginning of batch-reaction calculations.
-----------------------------------------

Reaction step 1.

Using solution 1. Tank2 Seawater
Using reaction 1.

Reaction 1.

  2.322e+00 moles of the following reaction have been added:

                 Relative
Reactant            moles

NaOH                 1.00000

                 Relative
Element             moles
H                    1.00000
Na                   1.00000
O                    1.00000

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

C                 2.354e-03   1.489e+01
Ca                1.066e-02   6.746e+01
Cl                5.657e-01   3.580e+03
Fe                3.711e-08   2.349e-04
K                 1.058e-02   6.694e+01
Mg                5.507e-02   3.485e+02
Mn                3.773e-09   2.387e-05
N                 6.571e-06   4.158e-02
Na                4.858e-01   3.074e+03
S                 2.926e-02   1.851e+02
Si                7.382e-05   4.672e-01

----------------------------Description of solution----------------------------

                                       pH  =   8.585      Charge balance
                                       pe  =  12.528      Adjusted to redox equilibrium
      Specific Conductance (µS/cm,  19°C)  = 46699
                          Density (g/cm³)  =   1.02484
                               Volume (L)  = 6399.24350
                        Activity of water  =   0.981
                 Ionic strength (mol/kgw)  =   6.769e-01
                       Mass of water (kg)  =   6.328e+03
                 Total alkalinity (eq/kg)  =   2.772e-03
                       Total CO2 (mol/kg)  =   2.354e-03
                         Temperature (°C)  =  19.34
                  Electrical balance (eq)  =   5.009e+00
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =   0.07
                               Iterations  =  19
                                  Total H  = 7.025051e+05
                                  Total O  = 3.520357e+05
------------------
End of simulation.
------------------

-----------------------------------------
Beginning of batch-reaction calculations.
-----------------------------------------

Reaction step 1.

Using solution 2. Solution after simulation 2.
Using pure phase assemblage 1.

-------------------------------Phase assemblage--------------------------------

                                                      Moles in assemblage
Phase               SI  log IAP  log K(T, P)   Initial       Final       Delta

CO2(g)           -3.50    -4.90     -1.40    1.000e+01   9.063e+00  -9.374e-01

-----------------------------Solution composition------------------------------

Elements           Molality       Moles

C                 2.502e-03   1.583e+01
Ca                1.066e-02   6.746e+01
Cl                5.657e-01   3.580e+03
Fe                3.711e-08   2.349e-04
K                 1.058e-02   6.694e+01
Mg                5.507e-02   3.485e+02
Mn                3.773e-09   2.387e-05
N                 6.571e-06   4.158e-02
Na                4.858e-01   3.074e+03
S                 2.926e-02   1.851e+02
Si                7.382e-05   4.672e-01

----------------------------Description of solution----------------------------

                                       pH  =   8.352      Charge balance
                                       pe  =  12.761      Adjusted to redox equilibrium
      Specific Conductance (µS/cm,  19°C)  = 46702
                          Density (g/cm³)  =   1.02484
                               Volume (L)  = 6399.27095
                        Activity of water  =   0.981
                 Ionic strength (mol/kgw)  =   6.770e-01
                       Mass of water (kg)  =   6.328e+03
                 Total alkalinity (eq/kg)  =   2.772e-03
                       Total CO2 (mol/kg)  =   2.502e-03
                         Temperature (°C)  =  19.34
                  Electrical balance (eq)  =   5.009e+00
 Percent error, 100*(Cat-|An|)/(Cat+|An|)  =   0.07
                               Iterations  =   8
                                  Total H  = 7.025051e+05
                                  Total O  = 3.520376e+05

------------------
End of simulation.
------------------

It appears that in the first section, TA and DIC are recalculated as 2.408e-03 eq/kg and 2.354e-03 mol/kg, respectively. Where did the extra carbon come from, when my input was C(4) 2271.1 umol/kgs?

2) Any suggestions on the pe/ redox selections? How much does this matter with the phreeq database?

3) On the last step, I equilibrate with CO2 in the atmosphere. Is there any way to differentiate air temperature from solution temperature?

Thanks for your help! It's been really interesting delving into Phreeqc-- I wish I had encountered it earlier!

[dlparkhurst]

(1) PHREEQC units are usually mol/kgw; so there is a conversion from mol/kgs to mol/kgw. Basically, the mass of solutes is subtracted from a kilogram of solution to find the mass of water. The moles of solutes are then divided by the mass of water to obtain mol/kgw.

If you want to see mol/kgs, you can use USER_PRINT, USER_PUNCH, or USER_GRAPH and Basic functions to make the calculation (with phreeqc.dat, Amm.dat, or pitzer.dat databases). for example:

10 tdic_mol_kgs = TOTMOL("C(4)") / (RHO * SOLN_VOL)

(2) I'm not sure how you want to deal with ammonium. You could ignore it (leave it out of SOLUTION), you could use Amm.dat and define N(-3) as Amm, which would be treated as unreactive, or you could leave it N(-3), but it will be oxidized to N(5) by dissolved oxygen in the reaction calculation. Fe will be calculated to be Fe(3) and Mn will be Mn(2), which should result in negligible effects regardless of your choice of initial redox state.

(3) You state the CO2 concentration as 420 ppm, which would correspond to log10(PCO2) of -3.38 rather than -3.5.  I'm not sure how the air temperature would affect the calculation; presumably some boundary layer effects with a transition from air to water temperature. Without getting too complicated, I'd just go with the water temperature.

I have no idea how to convert to a seawater pH scale.