on serpentinization

[w1nterfell]

Hello all,

I am currently setting up an input file that should model serpentinization reactions and H2-gas production on saturns moon enceladus from 0 to 200°C and up to 100atm using the phreeqc.dat database. Before acutally allowing for any reactions to occur I am testing solubilities of gases at these Ts and Ps. I experimented with gas_phase and equilibrium_phases and am currently having a small issue here...

the setup is simple, a solution (the moons ocean) is first equilibrated with a gas_phase (its atmosphere). the resulting solution is then equilibrated with another gas phase allowing  and T is raised to 200°C. This appreas to work rather well, however when trying reaction_pressure on the model the pressure is not considered by phreeqc. Here is the input file:


phases 1
Fayalite
        Fe2SiO4 +4.0000 H+  =  + 1.0000 H4SiO4 + 2.0000 Fe++
        log_k           19.1113
   -delta_H   -152.256   kJ/mol   # Calculated enthalpy of reaction   Fayalite
   -vm 46.390 #supcrt92
#   Enthalpy of formation:   -354.119 kcal/mol
        -analytic 1.3853e+001 -3.5501e-003 7.1496e+003 -6.8710e+000 -6.3310e+004
#       -Range:  0-300
Forsterite
        Mg2SiO4 +4.0000 H+  =  + 1.0000 H4SiO4 + 2.0000 Mg++
        log_k           27.8626
   -delta_H   -205.614   kJ/mol   # Calculated enthalpy of reaction   Forsterite
   -vm 43.790 #supcrt92
#   Enthalpy of formation:   -520 kcal/mol
        -analytic -7.6195e+001 -1.4013e-002 1.4763e+004 2.5090e+001 -3.0379e+005
#       -Range:  0-300
Enstatite
        MgSiO3 +2.0000 H+ + 1.0000 H2O =   + 1.0000 Mg++ + 1.0000 H4SiO4
        log_k           11.3269
   -delta_H   -82.7302   kJ/mol   # Calculated enthalpy of reaction   Enstatite
   -vm 31.276 #supcrt92
#   Enthalpy of formation:   -369.686 kcal/mol
        -analytic -4.9278e+001 -3.2832e-003 9.5205e+003 1.4437e+001 -5.4324e+005
#       -Range:  0-300
Diopside
        CaMgSi2O6 +4.0000 H+ + 2.0000 H2O =  + 1.0000 Ca++ + 1.0000 Mg++  + 2.0000 H4SiO4
        log_k           20.9643
   -delta_H   -133.775   kJ/mol   # Calculated enthalpy of reaction   Diopside
   -vm 66.090 #supcrt92
#   Enthalpy of formation:   -765.378 kcal/mol
        -analytic 7.1240e+001 1.5514e-002 8.1437e+003 -3.0672e+001 -5.6880e+005
#       -Range:  0-300
Ferrosilite
        FeSiO3 +2.0000 H+ + 1.0000 H2O =  + 1.0000 Fe++  + 1.0000 H4SiO4
        log_k           7.4471
   -delta_H   -60.6011   kJ/mol   # Calculated enthalpy of reaction   Ferrosilite
   -vm 32.952 #supcrt92
#   Enthalpy of formation:   -285.658 kcal/mol
        -analytic 9.0041e+000 3.7917e-003 5.1625e+003 -6.3009e+000 -3.9565e+005
#       -Range:  0-300
Magnetite
        Fe3O4 +8.0000 H+  =  + 1.0000 Fe++ + 2.0000 Fe+++ + 4.0000 H2O
        log_k           10.4724
   -delta_H   -216.597   kJ/mol   # Calculated enthalpy of reaction   Magnetite
   -Vm 44.524 #supcrt92
   #   Enthalpy of formation:   -267.25 kcal/mol
        -analytic -3.0510e+002 -7.9919e-002 1.8709e+004 1.1178e+002 2.9203e+002
#       -Range:  0-300

Brucite
        Mg(OH)2 +2.0000 H+  =  + 1.0000 Mg++ + 2.0000 H2O
        log_k           16.2980
   -delta_H   -111.34   kJ/mol   # Calculated enthalpy of reaction   Brucite
   -vm 24.630 #supcrt92
   #   Enthalpy of formation:   -221.39 kcal/mol
        -analytic -1.0280e+002 -1.9759e-002 9.0180e+003 3.8282e+001 1.4075e+002
#       -Range:  0-300
Talc
        Mg3Si4O10(OH)2 + 4 H2O + 6 H+ = 3 Mg+2 + 4 H4SiO4
        log_k           21.1383
   -delta_H   -148.737   kJ/mol   # Calculated enthalpy of reaction   Talc
#   Enthalpy of formation:   -1410.92 kcal/mol
        -analytic 1.1164e+001 2.4724e-002 1.9810e+004 -1.7568e+001 -1.8241e+006
#       -Range:  0-300
Antigorite
#        Mg48Si24O85(OH)62 +96.0000 H+  =  + 34.0000 SiO2 + 48.0000 Mg++ + 79.0000 H2O
        Mg48Si34O85(OH)62 +96.0000 H+  =  + 34.0000 H4SiO4 + 48.0000 Mg++ + 11 H2O
        log_k           477.1943
   -delta_H   -3364.43   kJ/mol   # Calculated enthalpy of reaction   Antigorite
#   Enthalpy of formation:   -17070.9 kcal/mol
        -analytic -8.1630e+002 -6.7780e-002 2.5998e+005 2.2029e+002 -9.3275e+006

Prehnite
        Ca2Al2Si3O10(OH)2 +10.0000 H+  =  + 2.0000 Al+++ + 2.0000 Ca++ + 3.0000 H4SiO4
        log_k           32.9305
   -delta_H   -311.875   kJ/mol   # Calculated enthalpy of reaction   Prehnite
#   Enthalpy of formation:   -1481.65 kcal/mol
        -analytic -3.5763e+001 -2.1396e-002 2.0167e+004 6.3554e+000 -7.4967e+005
#       -Range:  0-300   
Greenalite
        Fe3Si2O5(OH)4 +6.0000 H+  =  + 2.0000 H4SiO4 + 3.0000 Fe++ + 1.0000 H2O
        log_k           22.6701
   -delta_H   -165.297   kJ/mol   # Calculated enthalpy of reaction   Greenalite
#   Enthalpy of formation:   -787.778 kcal/mol
        -analytic -1.4187e+001 -3.8377e-003 1.1710e+004 1.6442e+000 -4.8290e+005
#       -Range:  0-300
Chrysotile
        Mg3Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 3 Mg+2
        log_k           31.1254
   -delta_H   -218.041   kJ/mol   # Calculated enthalpy of reaction   Chrysotile
#   Enthalpy of formation:   -1043.12 kcal/mol
        -analytic -9.2462e+001 -1.1359e-002 1.8312e+004 2.9289e+001 -6.2342e+005
#       -Range:  0-300
Quartz
        SiO2 + 2 H2O = H4SiO4
        log_k           -3.9993
   -delta_H   32.949   kJ/mol   # Calculated enthalpy of reaction   Quartz
#   Enthalpy of formation:   -217.65 kcal/mol
        -analytic 7.7698e-002 1.0612e-002 3.4651e+003 -4.3551e+000 -7.2138e+005
#       -Range:  0-300

SOLUTION 1   #composition of surface ocean
 -temp 0
Mg 0.025
N(0) 0.011
Ca 0.01
Fe 1.51e-9
Si 1.61e-4
C(-4) 9.1e-3
C(+4) 9e-3
S(-2) 2.1e-5
gas_phase  1    #composition fo gas plume
-fixed_volume
-volume 1
-temp 0
H2O(g)   0.9
CO2(g)   0.053
h2(g)    0.39
n2(g)    0.011
ch4(g)   9.1e-3
H2S(g)   2.1e-5
-equilibrate 1
save solution 1
end
use solution 1
gas_phase 2
-fixed_volume
-volume 1
-temp 0
h2(g)   0
N2(g)   0
ch4(g)  0
H2S(g)  0
H2O(g)  0
CO2(g)  0
REACTION_temperature
0 200 in 9
USER_GRAPH 1
 -axis_titles "temperature / oC" "Gas production / moles"
 -axis_scale x_axis 0 200
 -headings Gases H2(g) N2(g) CH4(g) CO2(g)
 -connect_simulations false
 10 plot_xy -1, 0, symbol = None , line_width = 0
 20 plot_xy tc, gas("h2(g)"), color = Red, symbol = Circle
 30 plot_xy tc, gas("N2(g)"), color = Blue, symbol = Circle
 40 plot_xy tc, gas("ch4(g)"), color = green, symbol = Plus
 #50 plot_xy tc, gas("H2o(g)"), color = green, symbol = Plus
 60 plot_xy tc, gas("co2(g)"), color = black, symbol = square
 END

The goal  is to add rocks of different composition using REACTION to the solution in order to estimate H2-gas production at varying minrealogies, Ts and Ps, and compare these to results attained from a similar model setup that uses EQ3/6.
After reading through the manual gas_phases using a fixed volume should be the correct approach, but how to consider higher pressures? Any suggestions, hints and improvements on the concept are much appreciated :).

[dlparkhurst]

The pressure in the GAS_PHASE definition is only used to calculate the number of moles of gases in the gas phase. Roughly PV = nRT, where P, V, and T are specified in the GAS_PHASE definition.

The reaction pressure is by default defined by the pressure of the solution that is used in the reaction. If you want to set the pressure of the reaction, then you should use the REACTION_PRESSURE keyword.

[w1nterfell]

thanks for the quick reply.

I have in fact tried the exact same input file using reaction_pressure instead of reaction_temp to estimate changes in gas phase composition from 1 to 1000atms. However, it appears that this keyword is ignored by gas_phase and solution, and the partial pressures within the phase do not change. Am I confusing something here?

regards

[dlparkhurst]

Sorry, I was not considering your gas phase carefully.

But what do you want? You cannot simultaneously fix the pressure and the volume of the gas phase. You are using a fixed volume GAS_PHASE. Pressure will be calculated as a function of reactions and temperature.

If you want to impose a pressure on the gas phase, then it must be a fixed-pressure gas phase, rather than a fixed volume.

Perhaps you do not want a gas phase at all, but simply look at the calculated partial pressures as a function of reaction, temperature, and pressure.

[w1nterfell]

Ah...well the goal is to test H2-gas production at different pressures up to around 200atms. So as I understand it a fixed volume gas phase is used and pressure by reaction_pressure is increased then this pressure is exerted on the solution, which then calculates moles of gases that go into the gas phase. However, in my experiment, this does not happen. Neither does the pressure of the gas_phase change (you accounted for that, thank you :)), but the pressure of the solution does not change either. That is what is befuddles me...it seems as if the reaction_pressure keyword is ignored.

"Perhaps you do not want a gas phase at all,"

not quite, there should be a gas phase of a known composition of gases present at all times (an "atmosphere" above the solution into which gases can be emitted from the solution or are dissolved into solution in response to reaction, pressure, and temperature changes)

"but simply look at the calculated partial pressures as a function of reaction, temperature, and pressure. "

exactly. I am trying to estimate changes in gas production/dissolution into of from "gas_phase" by changing these three boundary conditions.

here is the input file using phreeqc.dat:


phases 1 #from llnl.dat rewritten for phreeqc.dat
Fayalite
        Fe2SiO4 +4.0000 H+  =  + 1.0000 H4SiO4 + 2.0000 Fe++
        log_k           19.1113
   -delta_H   -152.256   kJ/mol   # Calculated enthalpy of reaction   Fayalite
   -vm 46.390 #supcrt92
#   Enthalpy of formation:   -354.119 kcal/mol
        -analytic 1.3853e+001 -3.5501e-003 7.1496e+003 -6.8710e+000 -6.3310e+004
#       -Range:  0-300
Forsterite
        Mg2SiO4 +4.0000 H+  =  + 1.0000 H4SiO4 + 2.0000 Mg++
        log_k           27.8626
   -delta_H   -205.614   kJ/mol   # Calculated enthalpy of reaction   Forsterite
   -vm 43.790 #supcrt92
#   Enthalpy of formation:   -520 kcal/mol
        -analytic -7.6195e+001 -1.4013e-002 1.4763e+004 2.5090e+001 -3.0379e+005
#       -Range:  0-300
Enstatite
        MgSiO3 +2.0000 H+ + 1.0000 H2O =   + 1.0000 Mg++ + 1.0000 H4SiO4
        log_k           11.3269
   -delta_H   -82.7302   kJ/mol   # Calculated enthalpy of reaction   Enstatite
   -vm 31.276 #supcrt92
#   Enthalpy of formation:   -369.686 kcal/mol
        -analytic -4.9278e+001 -3.2832e-003 9.5205e+003 1.4437e+001 -5.4324e+005
#       -Range:  0-300
Diopside
        CaMgSi2O6 +4.0000 H+ + 2.0000 H2O =  + 1.0000 Ca++ + 1.0000 Mg++  + 2.0000 H4SiO4
        log_k           20.9643
   -delta_H   -133.775   kJ/mol   # Calculated enthalpy of reaction   Diopside
   -vm 66.090 #supcrt92
#   Enthalpy of formation:   -765.378 kcal/mol
        -analytic 7.1240e+001 1.5514e-002 8.1437e+003 -3.0672e+001 -5.6880e+005
#       -Range:  0-300
Ferrosilite
        FeSiO3 +2.0000 H+ + 1.0000 H2O =  + 1.0000 Fe++  + 1.0000 H4SiO4
        log_k           7.4471
   -delta_H   -60.6011   kJ/mol   # Calculated enthalpy of reaction   Ferrosilite
   -vm 32.952 #supcrt92
#   Enthalpy of formation:   -285.658 kcal/mol
        -analytic 9.0041e+000 3.7917e-003 5.1625e+003 -6.3009e+000 -3.9565e+005
#       -Range:  0-300
Magnetite
        Fe3O4 +8.0000 H+  =  + 1.0000 Fe++ + 2.0000 Fe+++ + 4.0000 H2O
        log_k           10.4724
   -delta_H   -216.597   kJ/mol   # Calculated enthalpy of reaction   Magnetite
   -Vm 44.524 #supcrt92
   #   Enthalpy of formation:   -267.25 kcal/mol
        -analytic -3.0510e+002 -7.9919e-002 1.8709e+004 1.1178e+002 2.9203e+002
#       -Range:  0-300

Brucite
        Mg(OH)2 +2.0000 H+  =  + 1.0000 Mg++ + 2.0000 H2O
        log_k           16.2980
   -delta_H   -111.34   kJ/mol   # Calculated enthalpy of reaction   Brucite
   -vm 24.630 #supcrt92
   #   Enthalpy of formation:   -221.39 kcal/mol
        -analytic -1.0280e+002 -1.9759e-002 9.0180e+003 3.8282e+001 1.4075e+002
#       -Range:  0-300
Talc
        Mg3Si4O10(OH)2 + 4 H2O + 6 H+ = 3 Mg+2 + 4 H4SiO4
        log_k           21.1383
   -delta_H   -148.737   kJ/mol   # Calculated enthalpy of reaction   Talc
#   Enthalpy of formation:   -1410.92 kcal/mol
        -analytic 1.1164e+001 2.4724e-002 1.9810e+004 -1.7568e+001 -1.8241e+006
#       -Range:  0-300
Antigorite
#        Mg48Si24O85(OH)62 +96.0000 H+  =  + 34.0000 SiO2 + 48.0000 Mg++ + 79.0000 H2O
        Mg48Si34O85(OH)62 +96.0000 H+  =  + 34.0000 H4SiO4 + 48.0000 Mg++ + 11 H2O
        log_k           477.1943
   -delta_H   -3364.43   kJ/mol   # Calculated enthalpy of reaction   Antigorite
#   Enthalpy of formation:   -17070.9 kcal/mol
        -analytic -8.1630e+002 -6.7780e-002 2.5998e+005 2.2029e+002 -9.3275e+006

Prehnite
        Ca2Al2Si3O10(OH)2 +10.0000 H+  =  + 2.0000 Al+++ + 2.0000 Ca++ + 3.0000 H4SiO4
        log_k           32.9305
   -delta_H   -311.875   kJ/mol   # Calculated enthalpy of reaction   Prehnite
#   Enthalpy of formation:   -1481.65 kcal/mol
        -analytic -3.5763e+001 -2.1396e-002 2.0167e+004 6.3554e+000 -7.4967e+005
#       -Range:  0-300   
Greenalite
        Fe3Si2O5(OH)4 +6.0000 H+  =  + 2.0000 H4SiO4 + 3.0000 Fe++ + 1.0000 H2O
        log_k           22.6701
   -delta_H   -165.297   kJ/mol   # Calculated enthalpy of reaction   Greenalite
#   Enthalpy of formation:   -787.778 kcal/mol
        -analytic -1.4187e+001 -3.8377e-003 1.1710e+004 1.6442e+000 -4.8290e+005
#       -Range:  0-300
Chrysotile
        Mg3Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 3 Mg+2
        log_k           31.1254
   -delta_H   -218.041   kJ/mol   # Calculated enthalpy of reaction   Chrysotile
#   Enthalpy of formation:   -1043.12 kcal/mol
        -analytic -9.2462e+001 -1.1359e-002 1.8312e+004 2.9289e+001 -6.2342e+005
#       -Range:  0-300
Quartz
        SiO2 + 2 H2O = H4SiO4
        log_k           -3.9993
   -delta_H   32.949   kJ/mol   # Calculated enthalpy of reaction   Quartz
#   Enthalpy of formation:   -217.65 kcal/mol
        -analytic 7.7698e-002 1.0612e-002 3.4651e+003 -4.3551e+000 -7.2138e+005
#       -Range:  0-300

SOLUTION 1 #composition of surface ocean
 -temp 0
Mg 0.025
N(0) 0.011
Ca 0.01
Fe 1.51e-9
Si 1.61e-4
C(-4) 9.1e-3
C(+4) 9e-3
S(-2) 2.1e-5

gas_phase  1    #composition of atmosphere
-fixed_volume
-volume 1
-temp 0
H2O(g)   0.9
CO2(g)   0.053
h2(g)    0.39
n2(g)    0.011
ch4(g)   9.1e-3
H2S(g)   2.1e-5

save solution 1
end
use solution 1

gas_phase 2
-fixed_volume
-volume 1
-temp 0
H2O(g)   0.9
CO2(g)   0.053
h2(g)    0.39
n2(g)    0.011
ch4(g)   9.1e-3
H2S(g)   2.1e-5

reaction_pressure 1 # increasing pressure
1 1000 in 9

[dlparkhurst]

You cannot set the pressure for a fixed volume GAS_PHASE. That is just not the way it works. As you saw, the REACTION_PRESSURE is ignored and the pressure is determined by reactions between the gas phase and the solution.

But I still don't understand the physical situation. If the solution is in contact
with the atmosphere, it seems the pressure should remain equal to the atmospheric
pressure. Surely the atmosphere does not range from 1 to 100 (or 1000) atm. If the
pressure is generated with depth below moon surface then the appropriate gas phase
would be a fixed-pressure gas phase that would correspond to the pressure at a given
depth, and a gas phase would not exist unless the sum of the partial pressures of the
components were greater than the total pressure.

So walk me through this. The pressure is 1 atm; temperature 0. No additional H2 has
been created. What is the gas phase composition that you want? You have the
following:

GAS_PHASE 2
-fixed_pressure
-volume 1
-temp 0
H2O(g)   0.9
CO2(g)   0.053
H2(g)    0.39
N2(g)    0.011
CH4(g)   9.1e-3
H2S(g)   2.1e-5
END

However, the log partial pressures calculated for your SOLUTION 1 are all less than
-2. So, nearly all of your gas will dissolve into the water and the corresponding gas
phase pressure will be much less than 1. Note that the log partial pressure of H2O(g)
in equilibrium will not differ much from -2 because activities of water do not vary too far from 1.

Do you want the water to be in equilibrium with this gas phase (excepting the water
partial pressure)?

If so, you can either (1) make the gas phase with a very large
volume, or (2) fix the gas partial pressures with EQUILIBRIUM_PHASES and a large number of moles of each gas in EQUILIBRIUM_PHASES. For (1), the
total pressure will then be the sum of the partial pressures of the gas phase, which
will be less than 1 atm for the given definition. For (2), you can set the total pressure
to 1 atm with REACTION_PRESSURE.

If not, what do you want?

 

[w1nterfell]

Dead David,

thank you for your quick reply and explanations...

"So walk me through this."

I had in mind the concept that, as mentioned previously, a surface ocean is in contact with an atmosphere. the former composition is estimated, the second is composed of available literature data measured from a gas plume (which is why I focused on a fixed_volume gas phase enabling the presence of this atmosphere). In the first step I assume that they are not in equilibrium so i equilibrate them and then save them. Subsequently, serpentinization sets in using REACTION and EQUILIBRIUM_PHASES at varying Ts and Ps. I use the saved solution and, since the atmosphere is still there, the saved gas_phase with a fixed volume. I then expected to observe changes of gas composition within gas_phase as REACTION sets in at different Ts and Ps. However as you point out:

"If the pressure is generated with depth below moon surface then the appropriate gas phase would be a fixed-pressure gas phase that would correspond to the pressure at a given depth, and a gas phase would not exist unless the sum of the partial pressures of the components were greater than the total pressure."

this clairfies a lot...  there doesnt even have to be a gas phase..hmm. this poses somewhat of a problem since that would suggest that there shouldnt be any "atmosphere" or gas plume in the first place using my input data.

"Do you want the water to be in equilibrium with this gas phase (excepting the water partial pressure)? "

yes, because a gas phase should be present into which newly formed H2gas is emitted via serpentinization (using REACTION). I'll try using a fixed_pressure large-vol. gas phase and equilibrium_phases as you suggested and see what happens...

please excuse if what i am writing is somewhat confusing...i havent used gas_phase at all so far and although i have read the manual several times on this subject i am struggling to find the correct concept for this specific model.

many thanks and regards...

[dlparkhurst]

Before you get too far, there are other issues.

First, there is the problem with the partial pressure of water in the atmosphere. That would require a temperature close to 100 C to produce.

Second, the atmosphere is not in redox equilibrium. H2 and CO2 at those partial pressures would not be in equilibrium (much like O2 and N2 are not at thermodynamic equilibrium). So, for starters, I suggest that you decouple some of the redox reactions. Uncouple C(4) with C(-4) by setting a much lower log K in the CH4 SOLUTION_SPECIES definition. With the Amm.dat database, use Ntg for N2, and Mtg for methane.

Here is a first shot at a file, with some H2 production due to oxidation of iron metal. You can consider other options for the reactions and redox coupling.

SOLUTION_SPECIES
CO3-2 + 10H+ + 8e- = CH4 + 3H2O
    log_k     -100
    -dw       1.85e-09
END
SOLUTION 1
    temp      0
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    C(4)      0.009    
    Ca        0.01
    Fe        1.51e-09
    Mg        0.025
    Mtg       0.0091   
    Ntg       0.011   
    S(-2)     2.1e-05   
    Si        0.000161   
    H(0)      1      
    -water    1 # kg
END
EQUILIBRIUM_PHASES 1
    CO2(g)    -1.27572413  10
    H2(g)     -0.408935393 10
    H2S(g)    -4.677780705 10
    Mtg(g)    -2.040958608 10
    Ntg(g)    -1.958607315 10
END
USE solution 1
USE equilibrium_phases 1
SAVE solution 1
END
USE solution 1
USE equilibrium_phases 1
REACTION 1
    Fe         1
    1 millimoles
END

[w1nterfell]

wow, thank you for your suggestions and the input file :).

at this point i do have a question though..how is it that you can simply "decouple some of the redox reactions." in phreeqc without influencing (in this case) the carbon system or - for that matter - any other system? I mean if you simply change the log_k of any reaction doesnt that affect the distribution of all other species? So you chose to uncouple the lowest and highest redox state of carbon as CO2(g) is not in equilibrium...but you could also arbitrarily decouple f.e. S(6) and S(-2)?. I dont quite understand the logic behind the uncoupling of redox reactions in phreeqc.

In this regard I also dont quite understand the point of these "redox-uncoupled gases"...are mtg and ntg independent of solution composition, pH and pe, whereas H2(g) is not?

many thanks and regards

[dlparkhurst]

Redox equilibrium applies among all elements that have multiple redox states.

For your case, let's consider redox in the example that I sent.

Initially, for C there are C(4) and C(-4) redox states defined in Amm.dat. To decouple a redox state, say C(-4), it is necessary to remove or disable all of the species of that redox state defined in the database (Amm.dat). For C(-4) it is relatively easy because there is only one aqueous species and one mineral (CH4(g)). You can remove aqueous CH4 from Amm.dat, or you can make its concentration negligibly small by changing the K. You then need to redefine the C(-4) species as another "element". Fortunately, Mtg is already defined, so we can simply use that. The result is that C effectively has only one redox state, so it is independent of any redox reactions. If a redox state has a lot of species and minerals, it can be a lot of work. Fe(2) is decoupled from Fe(3) in one of the examples in the manual.

In my example, no N was defined in the solution (or minerals), so none of those redox reactions come into play. Zero-valent nitrogen was defined as Ntg.

For S, redox states and aqueous species are defined for S(-2) and S(6).

Fe has redox states and aqueous species for Fe(2) and Fe(3).

O(0) and H(0) are also defined in addition to water (H(1) and O(-2)).

One way to look at redox is that you can compute a pe from any two redox states of the same element. For redox equilibrium, the pe that you calculate from any pair will be the same as any other pair, so that pe(H(1)/H(0)) = pe(O(0)/O(-2)) = pe(Fe(2)/Fe(3)) = Pe(S(6)/S(-2)). This means that for any reaction, H, O, Fe, and S will adjust among their redox states to produce this equal pe.

Meanwhile C and Mtg will not react because they are not redox states of the same "element" (C and Mtg are as different as Ca and Mg). Any conversion between these two entities would have to be done with REACTION or KINETICS, for example

REACTION
CH4 -1
Ntg +1
1 umol

would convert some C to Ntg.

The reason to decouple C(4) from C(-4) is that by redox equilbrium, I think you would have no C(4) species. Try it. Leave the K for CH4 unchanged and add methane as C(-4). You then requilre pe(C(4)/C(-4)) to be equal to all the others, and carbon can be shuffled between the two redox states to obtain equilibrium.

[w1nterfell]

Dear David,

Thank you for your explanations. I have experimented with this setup for a while now and have left the SOLUTION_SPECIES as you defined it in the input file you posted above. I have also considered modeling this with and without a gas phase, and using pH-dependent dissolution kinetics of olivine, enstatite and diopside together with TRANSPORT to get a H2-production rate estimate. However, I would like to ask why you added "H(0) 1" to the input solution? Ideally there should be no (or as good as no) H2 in solution to begin with, and only be produced in the course of olivine dissolution and Fe oxidation. The question being if the H2 measured in the moon plume could have originated from such reactions or that it has a different origin independent of serpentinization.

many thanks and regards

[dlparkhurst]

I suspect it does not matter too much what you pick for H(0) in solution given that you are equilibrating with H2(g) in EQUILIBRIUM_PHASES. EQUILIBRIUM_PHASES will then be a source/sink for H2. In addition, the iron metal will generate H2. But to be honest, I don't know why I set H(0) in the solution other than it is a rough guess as to what the H2(aq) concentration will be when equilibrated with H2(g) in EQUILIBRIUM_PHASES.

[cchristophi]

What is the source of the data for the phases ?

[dlparkhurst]

If you are asking about thermodynamic data for PHASES in phreeqc.dat, it is mostly from WATEQ4F. Sources of data are listed in the manual for WATEQ4F, http://wwwbrr.cr.usgs.gov/projects/GWC_chemtherm/software.htm.

[cchristophi]

Yes thank you, I'm aware of that but I can't find all of them in there, Antigorite for example

[dlparkhurst]

I don't think wateq4f.dat or phreeqc.dat include antigorite.

We simply translated the data for llnl.dat from the Lawrence Livermore data set used in Geochemist's Workbench (thermo.com.V8.R6.230). You will have to try to find sources either from LLNL or Geochemist's Workbench publications.

The sit.dat database provides a reference for antigonite.

[cchristophi]

Thank you very much

[cchristophi]

Hello everyone. Does any of you have or know where I can get thermodynamic data on Bytownite?