Reactive transport with biogeochemistry: aquifer zoning

Description of biogeochemistry

The two species of microbe are: sulfate reducers that catalyse the reaction (1) and methanogens that catalyse the reaction (2) Acetate (CHCOO) is added to the aquifer from surrounding rocks at a rate of mol.m.yr (the m in the denominator is a volume of rock, not a volume of groundwater) which is a rate consistent with observations. This feeds the microbes. The initial microbial population is small.

The standard MOOSE geochemistry database contains the necessary entries (see the biogeochemistry page for more discussion):

    "methanogen": {
      "species": {
        "CH3COO-": -1.0,
        "H2O": -1.0,
        "CH4(aq)": 1.0,
        "HCO3-": 1.0
      },
      "molar volume": 1,
      "molecular weight": 1E9,
      "logk": [2.727, 2.641, 2.699, 2.933, 3.28, 3.788, 4.336, 4.789]
    },
    "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]
    },

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

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

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

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

• and (units: mol.kg) are the molality of acetate and SO, respectively (these are free concentrations, not the bulk compositions);

• mol.kg and mol.kg are the half-saturation constant of the acetate and SO, respectively;

• 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.

In addition, the growth yield of the sulfate reducer is assumed to be 4.3g.mol. The sulfate_reducer has a molar mass of 1000g/mol, so kinetic_biological_efficiency = 4.3E-3.

[UserObjects]
  [rate_sulfate_reducer]
    type = GeochemistryKineticRate
    kinetic_species_name = "sulfate_reducer"
    intrinsic_rate_constant = 31.536 # 1E-9 mol(acetate)/mg(biomass)/s = 31.536 mol(acetate)/g(biomass)/year
    multiply_by_mass = true
    promoting_species_names = 'CH3COO- SO4--'
    promoting_indices = '1 1'
    promoting_monod_indices = '1 1'
    promoting_half_saturation = '70E-6 200E-6'
    direction = dissolution
    kinetic_biological_efficiency = 4.3E-3 # 4.3 g(biomass)/mol(acetate) = 4.3E-3 mol(biomass)/mol(acetate) (because sulfate_reducer has molar mass of 1E3 g/mol)
    energy_captured = 45E3
    theta = 0.2
    eta = 1
  []
[]

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

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

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

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

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

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

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

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

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

• J.K.mol is the gas constant;

• (units: K) is the temperature.

In addition, the growth yield of the sulfate reducer is assumed to be 2g.mol. The methanogen has a molar mass of g/mol, so kinetic_biological_efficiency = 2.0E-9.

[UserObjects]
  [rate_methanogen]
    type = GeochemistryKineticRate
    kinetic_species_name = "methanogen"
    intrinsic_rate_constant = 63.072 # 2E-9 mol(acetate)/mg(biomass)/s = 63.072 mol(acetate)/g(biomass)/year
    multiply_by_mass = true
    promoting_species_names = 'CH3COO-'
    promoting_indices = '1'
    promoting_monod_indices = '1'
    promoting_half_saturation = '20E-3'
    direction = dissolution
    kinetic_biological_efficiency = 2.0E-9 # 2 g(biomass)/mol(acetate) = 2E-9 mol(biomass)/mol(acetate)  (because methanogen has molar mass of 1E9 g/mol)
    energy_captured = 24E3
    theta = 0.5
    eta = 1
  []
[]

In addition, both species of microbe are assumed to die at a rate of s, which is implemented in the following GeochemistryKineticRate userobjects:

[UserObjects]
  [death_sulfate_reducer]
    type = GeochemistryKineticRate
    kinetic_species_name = "sulfate_reducer"
    intrinsic_rate_constant = 0.031536E-3 # 1E-9 g(biomass)/g(biomass)/s = 0.031536 g(biomass)/g(biomass)/year = 0.031536E-3 mol(biomass)/g(biomass)/year (because sulfate_reducer has molar mass of 1E3 g/mol)
    multiply_by_mass = true
    direction = death
    eta = 0.0
  []
[]

[UserObjects]
  [death_methanogen]
    type = GeochemistryKineticRate
    kinetic_species_name = "methanogen"
    intrinsic_rate_constant = 0.031536E-9 # 1E-9 g(biomass)/g(biomass)/s = 0.031536 g(biomass)/g(biomass)/year = 0.031536E-9 mol(biomass)/g(biomass)/year (because methanogen has molar mass of 1E9 g/mol)
    multiply_by_mass = true
    direction = death
    eta = 0.0
  []
[]

To avoid the accumulation of HS from the sulfate reducers, the aquifer is assumed to include small amounts of siderite minerals, which reaction to form mackinawite:

The remainder of the geochemistry is standard, except for the source_species_rates that encode the groundwater flow that are sent from the flow App described below:

rate_Ca_diffuse = 6.66667E-9 # 2E-6 mol.m^-3.yr^-1 = 2E-9 mol.litre^-1.yr^-1 divided by porosity of 0.3
rate_CH3COO_diffuse = 13.3333E-9 # 4E-6 mol.m^-3.yr^-1 = 4E-9 mol.litre^-1.yr^-1 divided by porosity of 0.3
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 10
    xmin = 0
    xmax = 200000
  []
[]

[GlobalParams]
  point = '100000 0 0'
  reactor = reactor
[]

[SpatialReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  swap_into_basis = 'Siderite'
  swap_out_of_basis = 'Fe++'
  prevent_precipitation = 'Pyrite Troilite'
  charge_balance_species = "HCO3-"
  constraint_species = "H2O              Ca++             HCO3-            SO4--            CH3COO-          HS-              CH4(aq)          Siderite         H+"
# ASSUME that 1 litre of solution initially contains:
  constraint_value = "  1.0              1E-3             2E-3             0.04E-3          1E-9             1E-9             1E-9             1               -7.5"
  constraint_meaning = "kg_solvent_water bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition free_mineral     log10activity"
  constraint_unit = "   kg               moles            moles            moles            moles            moles            moles            cm3            dimensionless"
  controlled_activity_name = 'H+'
  controlled_activity_value = 3.16227E-8 # this is pH=7.5
  kinetic_species_name = "sulfate_reducer methanogen"
  kinetic_species_initial_value = '1E-6 1E-6'
  kinetic_species_unit = 'mg mg'
  source_species_names = "H2O              Ca++                       SO4--            CH3COO-                        HS-              CH4(aq)       Fe++"
  source_species_rates = "rate_H2O_per_1l  rate_Ca_per_1l_with_source rate_SO4_per_1l  rate_CH3COO_per_1l_with_source rate_HS_per_1l rate_CH4_per_1l rate_Fe_per_1l"
  ramp_max_ionic_strength_initial = 1
  ramp_max_ionic_strength_subsequent = 1
  execute_console_output_on = ''
  solver_info = true
  evaluate_kinetic_rates_always = true
  adaptive_timestepping = true
  abs_tol = 1E-14
  precision = 16
[]

[UserObjects]
  [rate_sulfate_reducer]
    type = GeochemistryKineticRate
    kinetic_species_name = "sulfate_reducer"
    intrinsic_rate_constant = 31.536 # 1E-9 mol(acetate)/mg(biomass)/s = 31.536 mol(acetate)/g(biomass)/year
    multiply_by_mass = true
    promoting_species_names = 'CH3COO- SO4--'
    promoting_indices = '1 1'
    promoting_monod_indices = '1 1'
    promoting_half_saturation = '70E-6 200E-6'
    direction = dissolution
    kinetic_biological_efficiency = 4.3E-3 # 4.3 g(biomass)/mol(acetate) = 4.3E-3 mol(biomass)/mol(acetate) (because sulfate_reducer has molar mass of 1E3 g/mol)
    energy_captured = 45E3
    theta = 0.2
    eta = 1
  []
  [death_sulfate_reducer]
    type = GeochemistryKineticRate
    kinetic_species_name = "sulfate_reducer"
    intrinsic_rate_constant = 0.031536E-3 # 1E-9 g(biomass)/g(biomass)/s = 0.031536 g(biomass)/g(biomass)/year = 0.031536E-3 mol(biomass)/g(biomass)/year (because sulfate_reducer has molar mass of 1E3 g/mol)
    multiply_by_mass = true
    direction = death
    eta = 0.0
  []
  [rate_methanogen]
    type = GeochemistryKineticRate
    kinetic_species_name = "methanogen"
    intrinsic_rate_constant = 63.072 # 2E-9 mol(acetate)/mg(biomass)/s = 63.072 mol(acetate)/g(biomass)/year
    multiply_by_mass = true
    promoting_species_names = 'CH3COO-'
    promoting_indices = '1'
    promoting_monod_indices = '1'
    promoting_half_saturation = '20E-3'
    direction = dissolution
    kinetic_biological_efficiency = 2.0E-9 # 2 g(biomass)/mol(acetate) = 2E-9 mol(biomass)/mol(acetate)  (because methanogen has molar mass of 1E9 g/mol)
    energy_captured = 24E3
    theta = 0.5
    eta = 1
  []
  [death_methanogen]
    type = GeochemistryKineticRate
    kinetic_species_name = "methanogen"
    intrinsic_rate_constant = 0.031536E-9 # 1E-9 g(biomass)/g(biomass)/s = 0.031536 g(biomass)/g(biomass)/year = 0.031536E-9 mol(biomass)/g(biomass)/year (because methanogen has molar mass of 1E9 g/mol)
    multiply_by_mass = true
    direction = death
    eta = 0.0
  []
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O H+ CH3COO- CH4(aq) HS- Ca++ HCO3- SO4-- Fe++"
    kinetic_minerals = "sulfate_reducer methanogen"
    equilibrium_minerals = "*"
    kinetic_rate_descriptions = "rate_sulfate_reducer death_sulfate_reducer rate_methanogen death_methanogen"
  []
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = 'min(0.1 * (t + 1), 100)'
  []
  end_time = 20000
[]

[AuxVariables]
  [rate_H2O_per_1l] # change in H2O per 1 litre of aqueous solution that we consider at each node
  []
  [rate_CH3COO_per_1l] # change in CH3COO- per 1 litre of aqueous solution that we consider at each node
  []
  [rate_CH4_per_1l] # change in CH4(aq) per 1 litre of aqueous solution that we consider at each node
  []
  [rate_HS_per_1l] # change in HS- per 1 litre of aqueous solution that we consider at each node
  []
  [rate_Ca_per_1l] # change in Ca++ per 1 litre of aqueous solution that we consider at each node
  []
  [rate_SO4_per_1l] # change in SO4-- per 1 litre of aqueous solution that we consider at each node
  []
  [rate_Fe_per_1l] # change in Fe++ per 1 litre of aqueous solution that we consider at each node
  []

  [rate_CH3COO_per_1l_with_source] # change in CH3COO- per 1 litre of aqueous solution that we consider at each node, including the diffuse source
  []
  [rate_Ca_per_1l_with_source] # change in Ca per 1 litre of aqueous solution that we consider at each node, including the diffuse source
  []

  [transported_H2O]
  []
  [transported_CH3COO]
  []
  [transported_CH4]
  []
  [transported_HS]
  []
  [transported_Ca]
  []
  [transported_SO4]
  []
  [transported_Fe]
  []
[]

[AuxKernels]
  [rate_CH3COO_per_1l_with_source]
    type = ParsedAux
    args = 'rate_CH3COO_per_1l'
    variable = rate_CH3COO_per_1l_with_source
    function = 'rate_CH3COO_per_1l + ${rate_CH3COO_diffuse}'
    execute_on = 'timestep_begin timestep_end'
  []
  [rate_Ca_per_1l_with_source]
    type = ParsedAux
    args = 'rate_Ca_per_1l'
    variable = rate_Ca_per_1l_with_source
    function = 'rate_Ca_per_1l + ${rate_Ca_diffuse}'
    execute_on = 'timestep_begin timestep_end'
  []
  [transported_H2O]
    type = GeochemistryQuantityAux
    variable = transported_H2O
    species = H2O
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_end'
  []
  [transported_CH3COO]
    type = GeochemistryQuantityAux
    variable = transported_CH3COO
    species = "CH3COO-"
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_end'
  []
  [transported_CH4]
    type = GeochemistryQuantityAux
    variable = transported_CH4
    species = "CH4(aq)"
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_end'
  []
  [transported_HS]
    type = GeochemistryQuantityAux
    variable = transported_HS
    species = "HS-"
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_end'
  []
  [transported_Ca]
    type = GeochemistryQuantityAux
    variable = transported_Ca
    species = "Ca++"
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_end'
  []
  [transported_SO4]
    type = GeochemistryQuantityAux
    variable = transported_SO4
    species = "SO4--"
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_end'
  []
  [transported_Fe]
    type = GeochemistryQuantityAux
    variable = transported_Fe
    species = "Fe++"
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_end'
  []
[]

[Postprocessors]
  [time]
    type = TimePostprocessor
  []
[]

[VectorPostprocessors]
  [data]
    type = LineValueSampler
    start_point = '0 0 0'
    end_point = '200000 0 0'
    num_points = 501 # NOTE
    sort_by = x
    variable = 'transported_CH4 transported_CH3COO transported_SO4 free_mg_sulfate_reducer free_mg_methanogen'
  []
[]

[Outputs]
  exodus = true
  [csv]
    type = CSV
    time_step_interval = 10
    execute_on = 'INITIAL TIMESTEP_END FINAL'
  []
[]

Flow

The operator-split method is used to couple the flow with the biogeochemistry, as discussed in the theory page. Since porosity does not change in this simulation, and the groundwater velocity is specified, only the groundwater (not the surrounding rock) is considered. The flow is implemented using GeochemistryTimeDerivative and ConservativeAdvection Kernels, acting on the mole number per litre of groundwater of each of the chemical components that flow through the aquifer:

[Kernels]
  [dot_H2O]
    type = GeochemistryTimeDerivative
    variable = conc_H2O
    save_in = rate_H2O_times_vv
  []
  [dot_CH3COO]
    type = GeochemistryTimeDerivative
    variable = conc_CH3COO
    save_in = rate_CH3COO_times_vv
  []
  [dot_CH4]
    type = GeochemistryTimeDerivative
    variable = conc_CH4
    save_in = rate_CH4_times_vv
  []
  [dot_HS]
    type = GeochemistryTimeDerivative
    variable = conc_HS
    save_in = rate_HS_times_vv
  []
  [dot_Ca]
    type = GeochemistryTimeDerivative
    variable = conc_Ca
    save_in = rate_Ca_times_vv
  []
  [dot_SO4]
    type = GeochemistryTimeDerivative
    variable = conc_SO4
    save_in = rate_SO4_times_vv
  []
  [dot_Fe]
    type = GeochemistryTimeDerivative
    variable = conc_Fe
    save_in = rate_Fe_times_vv
  []
  [adv_H2O]
    type = ConservativeAdvection
    velocity = velocity
    upwinding_type = full
    variable = conc_H2O
  []
  [adv_CH3COO]
    type = ConservativeAdvection
    velocity = velocity
    upwinding_type = full
    variable = conc_CH3COO
  []
  [adv_CH4]
    type = ConservativeAdvection
    velocity = velocity
    upwinding_type = full
    variable = conc_CH4
  []
  [adv_HS]
    type = ConservativeAdvection
    velocity = velocity
    upwinding_type = full
    variable = conc_HS
  []
  [adv_Ca]
    type = ConservativeAdvection
    velocity = velocity
    upwinding_type = full
    variable = conc_Ca
  []
  [adv_SO4]
    type = ConservativeAdvection
    velocity = velocity
    upwinding_type = full
    variable = conc_SO4
  []
  [adv_Fe]
    type = ConservativeAdvection
    velocity = velocity
    upwinding_type = full
    variable = conc_Fe
  []
[]

The rate of change of each component is saved into an AuxVariable, for instance rate_H2O_times_vv. This has units mol.yr, because the Variables represent the number of moles per litre of groundwater of each component at each node. To convert it to mol.yr.litre(groundwater), a ParsedAux is used, such as:

[AuxKernels]
  [rate_H2O_auxk]
    type = ParsedAux
    variable = rate_H2O
    args = 'rate_H2O_times_vv nodal_void_volume'
    function = 'rate_H2O_times_vv / nodal_void_volume'
  []
[]

Fresh groundwater is added at the left, by boundary conditions such as:

[BCs]
  [inject_H2O]
    type = DirichletBC
    boundary = 'left right'
    variable = conc_H2O
    value = ${eqm_H2O}
  []
[]

At the end of each flow timestep, the geochemistry solver is run:

[MultiApps]
  [react]
    type = TransientMultiApp
    input_files = bio_zoning_conc.i
    clone_parent_mesh = true
    execute_on = 'timestep_end' # This is critical
  []
[]

with the flow rates from the flow solver acting as source terms to the geochemistry simulation:

[Transfers]
  [changes_due_to_flow]
    type = MultiAppCopyTransfer
    to_multi_app = react
    source_variable = 'rate_H2O rate_CH3COO rate_CH4 rate_HS rate_Ca rate_SO4 rate_Fe' # change in mole number at every node / dt
    variable = 'rate_H2O_per_1l rate_CH3COO_per_1l rate_CH4_per_1l rate_HS_per_1l rate_Ca_per_1l rate_SO4_per_1l rate_Fe_per_1l' # change in moles at every node / dt
  []
[]

After the geochemistry solver has completed its timestep, the new concentrations are retrieved from the geochemistry model in preparation for the next flow timestep:

[Transfers]
  [transported_moles_from_geochem]
    type = MultiAppCopyTransfer
    from_multi_app = react
    source_variable = 'transported_H2O transported_CH3COO transported_CH4 transported_HS transported_Ca transported_SO4 transported_Fe'
    variable = 'conc_H2O conc_CH3COO conc_CH4 conc_HS conc_Ca conc_SO4 conc_Fe'
  []
[]

Results

As seen in Figure 1, the simulation gradual evolves to steady-state in which the aquifer is partitioned into two distinct zones: the upstream zone is inhabited by the sulfate reducers, who drive the sulfate concentration to low values by around 100km from the fresh groundwater source (which sits at ). The downstream zone is inhabited by methanogens, who produce significant dissolved methane.

References