Energy available for microbial respiration

This example closely follows Section 7.4 of Bethke (2007) in which the energy available for microbial respiration is computed using the Nernst Eh values of various redox half-reactions. Bethke (2007) has several pages of interesting explanation that are not repeated here: we are focussing primarily on the MOOSE model and validating the MOOSE results.

The following assumptions are made.

  • The following alternate oxidation states are in redox disequilibrium:

    • CH3COO-

    • CH4(aq)

    • Fe+++

    • H2(aq)

    • HS-

    • NH4+

    • NO2-

  • The pH is 6

  • The fugacity of CO(g) is fixed at 0.01, and it is used in the basis in place of the aqueous species HCO

  • There is 1 free cm of the mineral Fe(OH)(ppd) and that it is used in the basis in place of the redox aqueous species Fe.

  • The bulk composition of the other basis species are:

    • SO4– 0.3963E-3 mol (mg.kg)

    • Ca++ 0.407E-3 mol (mg.kg)

    • Cl- 0.46E-3 mol (mg.kg this is the charge-balance species)

    • Na+ 0.473E-3 mol (mg.kg)

    • Mg++ 0.0895E-3 mol (mg.kg)

    • NO3- 3.508E-5 mol (mg.kg )

    • CH4(aq) 2.712E-5 mol (mg.kg)

    • NO2- 9.456E-6 mol (mg.kg)

    • CH3COO- 5.5526E-6 mol (mg.kg)

    • Fe++ 3.895E-6 mol (mg.kg)

    • NH4+ 6.029E-6 mol (mg.kg)

    • HS- 1.6445E-6 mol (mg.kg)

      - The free composition of the remaining basis species are:

    • O2(aq) 3.399E-6 free molal (mg.kg free)

    • H2(aq) 0.004E-6 free molal (mol.kg, free (not bulk))

  • All mineralisation reactions are ignored.

MOOSE input file

The MOOSE input file that corresponding to this model description is

[TimeIndependentReactionSolver<<<{"href": "../../../syntax/TimeIndependentReactionSolver/index.html"}>>>]
  model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object (you must create this UserObject yourself)"}>>> = definition
  temperature<<<{"description": "The temperature (degC) of the aqueous solution"}>>> = 25
  swap_out_of_basis<<<{"description": "Species that should be removed from the model_definition's basis and be replaced with the swap_into_basis species"}>>> = 'HCO3-  Fe+++'
  swap_into_basis<<<{"description": "Species that should be removed from the model_definition's equilibrium species list and added to the basis.  There must be the same number of species in swap_out_of_basis and swap_into_basis.  These swaps are performed before any other computations during the initial problem setup. If this list contains more than one species, the swapping is performed one-by-one, starting with the first pair (swap_out_of_basis[0] and swap_into_basis[0]), then the next pair, etc"}>>> = '  CO2(g) Fe(OH)3(ppd)'
  charge_balance_species<<<{"description": "Charge balance will be enforced on this basis species.  This means that its bulk mole number may be changed from the initial value you provide in order to ensure charge neutrality.  After the initial swaps have been performed, this must be in the basis, and it must be provided with a bulk_composition constraint_meaning."}>>> = "Cl-"
# TDS = 80.5mg/g (roughly) = (mass_non-water) / (mass_solvent_water + mass_non-water),
# so with mass_solvent_water = 1kg, mass_non-water = 87.6mg and total_mass = 1.0876
  constraint_species<<<{"description": "Names of the species that have their values fixed to constraint_value with meaning constraint_meaning.  All basis species (after swap_into_basis and swap_out_of_basis) must be provided with exactly one constraint.  These constraints are used to compute the configuration during the initial problem setup, and in time-dependent simulations they may be modified as time progresses."}>>> = "H2O              H+          CO2(g)        Fe(OH)3(ppd) Cl-              Na+              Ca++             Mg++             SO4--            Fe++             H2(aq)             HS-              O2(aq)             CH4(aq)           NO3-              NO2-              NH4+              CH3COO-"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              1E-6        0.01          0.02914      0.46E-3          0.473E-3         0.407E-3         0.0895E-3        0.3963E-3        3.895E-6         0.004E-6           1.6445E-6        3.399E-6           2.712E-5          3.508E-5          9.456E-6          6.029E-6          5.5526E-6"
  constraint_meaning<<<{"description": "Meanings of the numerical values given in constraint_value.  kg_solvent_water: can only be applied to H2O and units must be kg.  bulk_composition: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species but not in kinetic species, and units must be moles or mass (kg, g, etc).  bulk_composition_with_kinetic: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species and in kinetic species, and units must be moles or mass (kg, g, etc).  free_concentration: can be applied to all basis species that are not gas and not H2O and not mineral, and represents the total amount of the basis species existing freely (not as secondary species) within the solution, and units must be molal or mass_per_kg_solvent.  free_mineral: can be applied to all mineral basis species, and represents the total amount of the mineral existing freely (precipitated) within the solution, and units must be moles, mass or cm3.  activity and log10activity: can be applied to basis species that are not gas and not mineral and not sorbing sites, and represents the activity of the basis species (recall pH = -log10activity), and units must be dimensionless.  fugacity and log10fugacity: can be applied to gases, and units must be dimensionless"}>>> = "kg_solvent_water activity    fugacity      free_mineral bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition free_concentration bulk_composition free_concentration bulk_composition  bulk_composition  bulk_composition  bulk_composition  bulk_composition"
  constraint_unit<<<{"description": "Units of the numerical values given in constraint_value.  Dimensionless: should only be used for activity or fugacity constraints.  Moles: mole number.  Molal: moles per kg solvent water.  kg: kilograms.  g: grams.  mg: milligrams.  ug: micrograms.  kg_per_kg_solvent: kilograms per kg solvent water.  g_per_kg_solvent: grams per kg solvent water.  mg_per_kg_solvent: milligrams per kg solvent water.  ug_per_kg_solvent: micrograms per kg solvent water.  cm3: cubic centimeters"}>>> = "   kg             dimensionless dimensionless moles        moles            moles            moles            moles            moles            moles            molal              moles            molal              moles             moles             moles             moles             moles"
  ramp_max_ionic_strength_initial<<<{"description": "The number of iterations over which to progressively increase the maximum ionic strength (from zero to max_ionic_strength) during the initial equilibration.  Increasing this can help in convergence of the Newton process, at the cost of spending more time finding the aqueous configuration."}>>> = 0 # not needed in this simple problem
  max_initial_residual<<<{"description": "Attempt to alter the initial-guess molalities so that the initial residual for the Newton process (that finds the aqueous configuration) is less than this number of moles"}>>> = 1E-2
  stoichiometric_ionic_str_using_Cl_only<<<{"description": "If set to true, the stoichiometric ionic strength will be set equal to Cl- molality (or max_ionic_strength if the Cl- molality is too big).  This flag overrides ionic_str_using_basis_molality_only"}>>> = true # for comparison with GWB
  mol_cutoff<<<{"description": "Information regarding species with molalities less than this amount will not be outputted"}>>> = 1E-5
  abs_tol<<<{"description": "If the residual of the algebraic system (measured in mol) is lower than this value, the Newton process (that finds the aqueous configuration) is deemed to have converged"}>>> = 1E-17
[]

[UserObjects<<<{"href": "../../../syntax/UserObjects/index.html"}>>>]
  [definition]
    type = GeochemicalModelDefinition<<<{"description": "User object that parses a geochemical database file, and only retains information relevant to the current geochemical model", "href": "../../../source/userobjects/GeochemicalModelDefinition.html"}>>>
    database_file<<<{"description": "The name of the geochemical database file"}>>> = "../../../database/moose_geochemdb.json"
    basis_species<<<{"description": "A list of basis components relevant to the aqueous-equilibrium problem. H2O must appear first in this list.  These components must be chosen from the 'basis species' in the database, the sorbing sites (if any) and the decoupled redox states that are in disequilibrium (if any)."}>>> = "H2O H+ Cl- Na+ Ca++ Mg++ SO4-- Fe++ H2(aq) HS- O2(aq) Fe+++ HCO3- CH3COO- CH4(aq) NH4+ NO2- NO3-"
    equilibrium_minerals<<<{"description": "A list of minerals that are in equilibrium with the aqueous solution.  All members of this list must be in the 'minerals' section of the database file"}>>> = "Fe(OH)3(ppd)"
    equilibrium_gases<<<{"description": "A list of gases that are in equilibrium with the aqueous solution and can have their fugacities fixed, at least for some time and spatial location.  All members of this list must be in the 'gas' section of the database file"}>>> = "CO2(g)"
    piecewise_linear_interpolation<<<{"description": "If true then use a piecewise-linear interpolation of logK and Debye-Huckel parameters, regardless of the interpolation type specified in the database file.  This can be useful for comparing with results using other geochemistry software"}>>> = true # for comparison with GWB
  []
[]
(modules/geochemistry/test/tests/redox_disequilibrium/microbial.i)

GWB

The Geochemists Workbench input file corresponding to this model is

# React script that is equivalent to the microbial.i MOOSE input file
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 25 C
decouple CH3COO-
decouple CH4(aq)
decouple Fe+++
decouple H2(aq)
decouple HS-
decouple NH4+
decouple NO2-
H2O          = 1 free kg
Cl-          = 0.46E-3 mol
balance on Cl-
H+           = 6 pH
swap CO2(g) for HCO3-
CO2(g)       = 0.01 fugacity
swap Fe(OH)3(ppd) for Fe+++
Fe(OH)3(ppd) = 0.02914 free mol
Na+          = 0.473E-3 mol
Ca++         = 0.407E-3 mol
Mg++         = 0.0895E-3 mol
SO4--        = 0.3963E-3 mol
Fe++         = 3.895E-6 mol
H2(aq)       = 0.004E-6 free molal
HS-          = 1.6445E-6 mol
O2(aq)       = 3.399E-6 free molal
CH4(aq)      = 2.712E-5 mol
NO3-         = 3.508E-5 mol
NO2-         = 9.456E-6 mol
NH4+         = 6.029E-6 mol
CH3COO-      = 5.5526E-6 mol
printout  species = long
suppress all
unsuppress Fe(OH)3(ppd)
epsilon = 1e-14
go
(modules/geochemistry/test/tests/redox_disequilibrium/microbial.rea)

Nernst potentials

Bethke (2007) computes the Nernst potentials for each redox couple as shown in Table 1

Table 1: Nernst potentials for redox couple

ReactionEh (mV)
e + O(aq) + H HO836
2e + 2H + NO HO + NO481
8e + 10H + NO 3HO + NH443
6e + 8H + NO 2HO + NH430
e + Fe Fe321
8e + 9H + SO 4HO + HS-126
4e + H + CHCOO HO + CH(aq)-145
8e + 9H + HCO 3HO + CH(aq)-188
2e + 2H H(aq)-199
8e + 9H + 2HCO 4HO + CHCOO-230

This output is produced by the Geochemists Workbench software. The output produced by geochemistry looks a little different:


Nernst potentials:
e- = 0.5*H2O - 1*H+ - 0.25*O2(aq);  Eh = 0.8361V
e- = 0.5*H2O - 1*H+ + 0.5*NO2- - 0.5*NO3-;  Eh = 0.4809V
e- = 0.375*H2O - 1.25*H+ + 0.125*NH4+ - 0.125*NO3-;  Eh = 0.4425V
e- = 1*Fe++ - 1*Fe+++;  Eh = 0.3206V
e- = 0.5*H2O - 1.125*H+ - 0.125*SO4-- + 0.125*HS-;  Eh = -0.1263V
e- = 0.375*H2O - 1.125*H+ - 0.125*HCO3- + 0.125*CH4(aq);  Eh = -0.1875V
e- = -1*H+ + 0.5*H2(aq);  Eh = -0.1986V
e- = 0.5*H2O - 1.125*H+ - 0.25*HCO3- + 0.125*CH3COO-;  Eh = -0.2298V

The differences are:

  • the reactions are all written with "e-" on the left-hand side

  • the values for Eh sometimes differ in their fourth significant figure (this is presumably due to geochemistry using higher precision than GWB)

  • the following reactions are missing:

    • 6e + 8H + NO 2HO + NH (Eh = 430) and

    • 4e + H + CHCOO HO + CH(aq) (Eh = -145)

The reason these reactions do not appear is that geochemistry outputs a minimal set only, from which others can be derived. For instance, the difference of the two equations produces the equation which produces the result mV when divided by 3. (Remember that Eh values quoted by GWB are normalised to the number of electrons, but when manipulating equations in the way just described, the number of electrons must be accounted for.)

References

  1. Craig M. Bethke. Geochemical and Biogeochemical Reaction Modeling. Cambridge University Press, 2 edition, 2007. doi:10.1017/CBO9780511619670.[BibTeX]