Biogeochemistry with a sulfate reducer

This follows the example in Section 18.5 of Bethke (2007). Suppose that a sulfate-reducing microbe acts in the presence of acetate (CHCOO) to catalyse the reaction (1) Bethke details the initial composition of the aqueous solution for this problem. Most important for this example are the initial bulk composition of the acetate, which is moles, the initial bulk composition of sulfate, which is large, and the initial mass of the sulfate reducer, which is 0.1mg.

In this example, microbe mortality is unimportant. For each mole turnover of this reaction, Bethke estimates the microbe mass increases by g.

The rate of Eq. (1) is assumed to be of the form (2) where

  • (units: kg) is the mass of solvent water;

  • mol.mg.smol.g.day is the intrinsic rate constant, where the mass (grams) in the denominator is the mass of microbe, and the direction is forward (dissolution of acetate only, not the reverse);

  • (units: g.kg) is the microbe concentration, where the mass (kg) in the denominator is the mass of solvent water;

  • (units: mol.kg) is the molality of CHCOO (this is the free concentration, not the bulk composition);

  • mol.kg is the half-saturation constant of the acetate;

  • is the activity product of Eq. (1);

  • is the equilibrium constant of Eq. (1);

  • J.mol is an estimate of the energy captured by the microbe for each mole turnover of Eq. (1);

  • J.K.mol is the gas constant;

  • (units: K) is the temperature.

As discussed on the main biogeochemistry page, there are two ways of modelling this scenario:

  1. the microbe is treated as a new basis species, and the acetate (CHCOO) as a kinetic species that is not an equilibrium with the aqueous solution;

  2. the microbe is treated as a kinetic species, and all other species (acetate, etc) are at equilibrium

Input files and results for both of these methods are found below.

Method 1

Here, the microbe is considered as a primary species. The aqueous solution as described by Bethke contains a number of ionic species: the important thing here is the existence of Biomass1 in the basis_species list:

[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 Na+ Ca++ Fe++ Cl- SO4-- HCO3- O2(aq) H+ Biomass1"
    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"}>>> = "Mackinawite" # other minerals make marginal difference
    kinetic_redox<<<{"description": "A list alternative oxidation states (eg Fe+++) whose dynamics are governed by a rate law.  These are not in equilibrium with the aqueous solution.  All members of this list must be in the 'redox couples' section of the database file."}>>> = "CH3COO-"
    kinetic_rate_descriptions<<<{"description": "A list of GeochemistryKineticRate UserObject names that define the kinetic rates.  If a kinetic species has no rate prescribed to it, its reaction rate will be zero"}>>> = "rate_sulfate_reducer"
    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 # comparison with GWB
  []
[]
(modules/geochemistry/test/tests/kinetics/bio_sulfate_1.i)

The Biomass1 species has a molar mass of g.mol.

The acetate is considered a kinetic species. Eq. (2) is implemented using the following GeochemistryKineticRate UserObject:

[UserObjects<<<{"href": "../../../syntax/UserObjects/index.html"}>>>]
  [rate_sulfate_reducer]
    type = GeochemistryKineticRate<<<{"description": "User object that defines a kinetic rate.  Note that more than one rate can be prescribed to a single kinetic_species: the sum the individual rates defines the overall rate.  GeochemistryKineticRate simply specifies the algebraic form for a kinetic rate: to actually use it in a calculation, you must use it in the GeochemicalModelDefinition.  The rate is intrinsic_rate_constant * area_quantity * (optionally, mass of kinetic_species in grams) * kinetic_molality^kinetic_molal_index / (kinetic_molality^kinetic_molal_index + kinetic_half_saturation^kinetic_molal_index)^kinetic_monod_index * (product_over_promoting_species m^promoting_index / (m^promoting_index + promoting_half_saturation^promiting_index)^promoting_monod_index) * |1 - (Q/K)^theta|^eta * exp(activation_energy / R * (1/T0 - 1/T)) * Direction(1 - (Q/K)).  Please see the markdown documentation for examples", "href": "../../../source/userobjects/GeochemistryKineticRate.html"}>>>
    kinetic_species_name<<<{"description": "The name of the kinetic species that will be controlled by this rate"}>>> = "CH3COO-"
    intrinsic_rate_constant<<<{"description": "The intrinsic rate constant for the reaction"}>>> = 0.0864 # 1E-9 mol/mg/s = 0.0864 mol/g/day
    multiply_by_mass<<<{"description": "Whether the rate should be multiplied by the kinetic_species mass (in grams)"}>>> = false
    kinetic_molal_index<<<{"description": "The rate is multiplied by kinetic_species_molality^kinetic_molal_index / (kinetic_species_molality^kinetic_molal_index + kinetic_half_saturation^kinetic_molal_index)^kinetic_monod_index"}>>> = 1.0
    kinetic_monod_index<<<{"description": "The rate is multiplied by kinetic_species_molality^kinetic_molal_index / (kinetic_species_molality^kinetic_molal_index + kinetic_half_saturation^kinetic_molal_index)^kinetic_monod_index"}>>> = 1.0
    kinetic_half_saturation<<<{"description": "The rate is multiplied by kinetic_species_molality^kinetic_molal_index / (kinetic_species_molality^kinetic_molal_index + kinetic_half_saturation^kinetic_molal_index)^kinetic_monod_index"}>>> = 70E-6
    promoting_species_names<<<{"description": "Names of any promoting species"}>>> = 'H2O Biomass1'
    promoting_indices<<<{"description": "Indices of the promoting species"}>>> = '1 1'
    direction<<<{"description": "Direction of reaction.  Let Q = the activity product of the kinetic reaction, and K = the equilibrium constant of the reaction.  Then direction means the following.  both = dissolution and precipitation are allowed.  (Specifically, if Q < K then dissolution will occur, that is, the kinetic species mass will decrease with time.  If Q > K then precipitation will occur, that is, the kinetic species mass will increase with time.)  dissolution = if Q < K then dissolution will occur, and when Q > K then the rate will be set to zero so that precipitation will be prevented.  precipitation = if Q > K then precipitation will occur, and when Q < K then the rate will be set to zero so that dissolution will be prevented.  raw = the rate will not depend on sgn(1 - (Q/K)), which means dissolution will occur if intrinsic_rate_constant > 0, and precipitation will occur when intrinsic_rate_constant < 0.  death = the rate will not depend on sgn(1 - (Q/K)), which means dissolution will occur if intrinsic_rate_constant > 0, and precipitation will occur when intrinsic_rate_constant < 0, and, in addition, no reactants will be produced or consumed by this kinetic reaction (only the kinetic species mass will change)."}>>> = dissolution
    non_kinetic_biological_catalyst<<<{"description": "Name of the primary or equilibrium species that acts as a biological catalyst."}>>> = Biomass1
    non_kinetic_biological_efficiency<<<{"description": "When one mole of the kinetic species dissolves, non_kinetic_biological_efficiency moles of the non_kinetic_biological_catalyst is created"}>>> = 4.3
    energy_captured<<<{"description": "In biologically-catalysed kinetic reactions, this is the energy captured by the cell, per mol of reaction turnover.  Specifically, for each mole of kinetic reaction, the microbe will produce m moles of ATP via a reaction such as ADP + PO4--- -> ATP + H2O, with free-energy change G (usually around 45 kJ/mol).  Then, energy_captured = m * G.  For non-biologically-catalysed reactions, this should be zero.  The impact of energy_captured is that the reaction's equilibrium constant is K_database * exp(-energy_captured / R / T_in_Kelvin)"}>>> = 45E3
    theta<<<{"description": "Theta parameter, which appears in |1 - (Q/K)^theta|^eta"}>>> = 0.2
    eta<<<{"description": "Eta parameter, which appears in |1 - (Q/K)^theta|^eta"}>>> = 1
  []
[]
(modules/geochemistry/test/tests/kinetics/bio_sulfate_1.i)

Most of the values are self-explanatory. The parameter direction = dissolution because the microbe acts to convert acetate into bicarbonate, not the reverse. Notice the non_kinetic_biological_catalyst. Because one mole of Eq. (1) produces of microbe,


non_kinetic_biological_efficiency = 4.3

Consider the initial composition. Above, it is specified that the initial bulk composition of acetate is moles. It is easy to use the geochemistry module to show that the molality of acetate is 0.0008643mol.kg (for instance, remove the sulfate_reducer from bio_sulfate_2.i mentioned below and use a TimeIndependentReactionSolver). This is the molality that enters into Eq. (2). There are two possibilities:

  • initialize the kinetic acetate to 0.0008643moles. Then the rate is likely to be reasonably accurate, but there won't be as much acetate to react compared with the case when moles are used, so the final biomass will be smaller.

  • initialize the kinetic acetate to moles. Then the rate will be slightly higher than it should be, but the full amount of acetate will react, so the final biomass will be reasonably correct. This approach is used to generate the graphs below.

This example illustrate the inaccuracies brought about by using the simple approach of treating the biomass as an equilibrium species that is generated by a kinetic reaction.

Method 2

Here, the microbe is considered as a kinetic species that is generated via Eq. (1). The other species (acetate, sulfate, etc) are considered primary species in equilibrium in the aqueous solution. The standard MOOSE geochemical database contains the required information:


    "sulfate_reducer": {
      "species": {
        "CH3COO-": -1.0,
        "SO4--": -1.0,
        "HCO3-": 2.0,
        "HS-": 1.0
      },
      "molar volume": 1,
      "molecular weight": 1E3,
	"logk": [8.502, 8.404, 8.451, 8.657, 9.023, 9.457, 9.917, 10.31]
    },

If this entry did not appear in the database, it could be manually entered, but the equilibrium constants would have to be derived from those already in the database (or from experimental measurements). This is easily done by using the GeochemicalModelInterrogator:

[GeochemicalModelInterrogator<<<{"href": "../../../syntax/GeochemicalModelInterrogator/index.html"}>>>]
  model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object"}>>> = definition
  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.  There must be the same number of species in swap_out_of_basis and swap_into_basis.  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"}>>> = 'O2(aq) HS-'
  swap_into_basis<<<{"description": "Species that should be removed from the model_definition's equilibrium species list and added to the basis"}>>> = 'HS- CH4(aq)'
  temperature<<<{"description": "Equilibrium constants will be computed at this temperature [degC].  This is ignored if interrogation=eqm_temperature."}>>> = 300.0
[]
(modules/geochemistry/test/tests/kinetics/bio_sulfate_0.i)

and the appropriate model definition:

[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 HCO3- SO4-- H+ O2(aq)"
    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
  []
[]
(modules/geochemistry/test/tests/kinetics/bio_sulfate_0.i)

After the entry exists in the database, the model description starts with defining the species present:

[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 Na+ Ca++ Fe++ Cl- SO4-- HCO3- O2(aq) H+ CH3COO-"
    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"}>>> = "Mackinawite" # other minerals make marginal difference
    kinetic_minerals<<<{"description": "A list of minerals whose dynamics are governed by a rate law.  These are not in equilibrium with the aqueous solution.  All members of this list must be in the 'minerals' section of the database file."}>>> = "sulfate_reducer"
    kinetic_rate_descriptions<<<{"description": "A list of GeochemistryKineticRate UserObject names that define the kinetic rates.  If a kinetic species has no rate prescribed to it, its reaction rate will be zero"}>>> = "rate_sulfate_reducer"
    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 # comparison with GWB
  []
[]
(modules/geochemistry/test/tests/kinetics/bio_sulfate_2.i)

The sulfate_reducer species has a molar mass of g.mol in the database file. It is unimportant that this is definitely incorrect: it is only important to remember this value when entering numerical values for the kinetic rate, below.

Eq. (2) is implemented using the following GeochemistryKineticRate UserObject:

[UserObjects<<<{"href": "../../../syntax/UserObjects/index.html"}>>>]
  [rate_sulfate_reducer]
    type = GeochemistryKineticRate<<<{"description": "User object that defines a kinetic rate.  Note that more than one rate can be prescribed to a single kinetic_species: the sum the individual rates defines the overall rate.  GeochemistryKineticRate simply specifies the algebraic form for a kinetic rate: to actually use it in a calculation, you must use it in the GeochemicalModelDefinition.  The rate is intrinsic_rate_constant * area_quantity * (optionally, mass of kinetic_species in grams) * kinetic_molality^kinetic_molal_index / (kinetic_molality^kinetic_molal_index + kinetic_half_saturation^kinetic_molal_index)^kinetic_monod_index * (product_over_promoting_species m^promoting_index / (m^promoting_index + promoting_half_saturation^promiting_index)^promoting_monod_index) * |1 - (Q/K)^theta|^eta * exp(activation_energy / R * (1/T0 - 1/T)) * Direction(1 - (Q/K)).  Please see the markdown documentation for examples", "href": "../../../source/userobjects/GeochemistryKineticRate.html"}>>>
    kinetic_species_name<<<{"description": "The name of the kinetic species that will be controlled by this rate"}>>> = "sulfate_reducer"
    intrinsic_rate_constant<<<{"description": "The intrinsic rate constant for the reaction"}>>> = 0.0864 # 1E-9 mol/mg/s = 0.0864 mol/g/day
    multiply_by_mass<<<{"description": "Whether the rate should be multiplied by the kinetic_species mass (in grams)"}>>> = true
    promoting_species_names<<<{"description": "Names of any promoting species"}>>> = 'CH3COO-'
    promoting_indices<<<{"description": "Indices of the promoting species"}>>> = 1
    promoting_monod_indices<<<{"description": "Indices of the monod denominators of the promoting species.  If not given, then the default is 0 for each promoting species, meaning that there is no monod form"}>>> = 1
    promoting_half_saturation<<<{"description": "Half-saturation constants for the monod expression.  If not given, then the default is 0 for each promoting species"}>>> = 70E-6
    direction<<<{"description": "Direction of reaction.  Let Q = the activity product of the kinetic reaction, and K = the equilibrium constant of the reaction.  Then direction means the following.  both = dissolution and precipitation are allowed.  (Specifically, if Q < K then dissolution will occur, that is, the kinetic species mass will decrease with time.  If Q > K then precipitation will occur, that is, the kinetic species mass will increase with time.)  dissolution = if Q < K then dissolution will occur, and when Q > K then the rate will be set to zero so that precipitation will be prevented.  precipitation = if Q > K then precipitation will occur, and when Q < K then the rate will be set to zero so that dissolution will be prevented.  raw = the rate will not depend on sgn(1 - (Q/K)), which means dissolution will occur if intrinsic_rate_constant > 0, and precipitation will occur when intrinsic_rate_constant < 0.  death = the rate will not depend on sgn(1 - (Q/K)), which means dissolution will occur if intrinsic_rate_constant > 0, and precipitation will occur when intrinsic_rate_constant < 0, and, in addition, no reactants will be produced or consumed by this kinetic reaction (only the kinetic species mass will change)."}>>> = both
    kinetic_biological_efficiency<<<{"description": "This is used when modelling biologically-catalysed reactions, when the biomass is treated as a kinetic species, and the reactants and reactant-products are in equilibrium in the aqueous solution.  When one mole of reaction is catalysed, the biomass increases by kinetic_biological_efficiency moles"}>>> = 4.3E-3
    energy_captured<<<{"description": "In biologically-catalysed kinetic reactions, this is the energy captured by the cell, per mol of reaction turnover.  Specifically, for each mole of kinetic reaction, the microbe will produce m moles of ATP via a reaction such as ADP + PO4--- -> ATP + H2O, with free-energy change G (usually around 45 kJ/mol).  Then, energy_captured = m * G.  For non-biologically-catalysed reactions, this should be zero.  The impact of energy_captured is that the reaction's equilibrium constant is K_database * exp(-energy_captured / R / T_in_Kelvin)"}>>> = 45E3
    theta<<<{"description": "Theta parameter, which appears in |1 - (Q/K)^theta|^eta"}>>> = 0.2
    eta<<<{"description": "Eta parameter, which appears in |1 - (Q/K)^theta|^eta"}>>> = 1
  []
[]
(modules/geochemistry/test/tests/kinetics/bio_sulfate_2.i)

Most of the values are self-explanatory. Note:

  • multiply_by_mass = true means the rate is multiplied by the mass (in grams) of the sulfate_reducer, as desired

  • there is no need to use H2O in the promoting_species because is just the mass of the sulfate_reducer in grams

  • kinetic_biological_efficiency = 4.3E-3 because this is the number of moles of sulfate_reducer that is created for one mole Eq. (1) turnover: recall that its molar mass is g.mol

Results

As expected, the two methods differ marginally: see Figure 1. Method 1 is slightly inaccurate because the full amount of acetate is used in the rate, instead of just the free molality as in method 2. Method 2 is initialized with the same bulk composition ( moles) but some of this is bound into secondary species, which slowly release the acetate into solution in order to react via Eq. (1).

Figure 1: Results of the biologically-catalysed sulfate reduction

References

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