Reactive transport with biogeochemistry: aquifer zoning
This follows the example in Section 33.2 of Bethke (2007). Groundwater flows through an aquifer from west to east (towards the positive "" direction). The aquifer is 200km long, has a porosity of 0.3, and the flow rate is 10m(water).m(rock).yr. The groundwater has pH 7.5, and contains 1mmolal Ca, 2mmolal HCO, 40molal SO, and small amounts of acetate, sulfide and methane. Living in the aquifer are two species of microbe: sulfate reducers and methanogens, which are discussed in detail below. The model definition is:
[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+ CH3COO- CH4(aq) HS- Ca++ HCO3- SO4-- Fe++"
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 methanogen"
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"}>>> = "*"
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 death_sulfate_reducer rate_methanogen death_methanogen"
[]
[]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<<<{"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"}>>> = 31.536 # 1E-9 mol(acetate)/mg(biomass)/s = 31.536 mol(acetate)/g(biomass)/year
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- SO4--'
promoting_indices<<<{"description": "Indices of the promoting species"}>>> = '1 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 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 200E-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)."}>>> = dissolution
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 # 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<<<{"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
[]
[]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<<<{"href": "../../../syntax/UserObjects/index.html"}>>>]
[rate_methanogen]
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"}>>> = "methanogen"
intrinsic_rate_constant<<<{"description": "The intrinsic rate constant for the reaction"}>>> = 63.072 # 2E-9 mol(acetate)/mg(biomass)/s = 63.072 mol(acetate)/g(biomass)/year
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"}>>> = '20E-3'
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
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"}>>> = 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<<<{"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)"}>>> = 24E3
theta<<<{"description": "Theta parameter, which appears in |1 - (Q/K)^theta|^eta"}>>> = 0.5
eta<<<{"description": "Eta parameter, which appears in |1 - (Q/K)^theta|^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<<<{"href": "../../../syntax/UserObjects/index.html"}>>>]
[death_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.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<<<{"description": "Whether the rate should be multiplied by the kinetic_species mass (in grams)"}>>> = true
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)."}>>> = death
eta<<<{"description": "Eta parameter, which appears in |1 - (Q/K)^theta|^eta"}>>> = 0.0
[]
[][UserObjects<<<{"href": "../../../syntax/UserObjects/index.html"}>>>]
[death_methanogen]
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"}>>> = "methanogen"
intrinsic_rate_constant<<<{"description": "The intrinsic rate constant for the reaction"}>>> = 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<<<{"description": "Whether the rate should be multiplied by the kinetic_species mass (in grams)"}>>> = true
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)."}>>> = death
eta<<<{"description": "Eta parameter, which appears in |1 - (Q/K)^theta|^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<<<{"href": "../../../syntax/Mesh/index.html"}>>>]
[gen]
type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 1
nx<<<{"description": "Number of elements in the X direction"}>>> = 10
xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = 0
xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 200000
[]
[]
[GlobalParams<<<{"href": "../../../syntax/GlobalParams/index.html"}>>>]
point = '100000 0 0'
reactor = reactor
[]
[SpatialReactionSolver<<<{"href": "../../../syntax/SpatialReactionSolver/index.html"}>>>]
model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object (you must create this UserObject yourself)"}>>> = definition
geochemistry_reactor_name<<<{"description": "The name that will be given to the GeochemistryReactor UserObject built by this action"}>>> = reactor
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"}>>> = 'Siderite'
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"}>>> = 'Fe++'
prevent_precipitation<<<{"description": "Mineral species in this list will be prevented from precipitating, irrespective of their saturation index (unless they are initially in the basis)"}>>> = 'Pyrite Troilite'
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."}>>> = "HCO3-"
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 Ca++ HCO3- SO4-- CH3COO- HS- CH4(aq) Siderite H+"
# ASSUME that 1 litre of solution initially contains:
constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = " 1.0 1E-3 2E-3 0.04E-3 1E-9 1E-9 1E-9 1 -7.5"
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 bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition free_mineral log10activity"
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 moles moles moles moles moles moles cm3 dimensionless"
controlled_activity_name<<<{"description": "The names of the species that have their activity or fugacity constrained. There should be an equal number of these names as values given in controlled_activity_value. NOTE: if these species are not in the basis, or they do not have an activity (or fugacity) constraint then their activity cannot be controlled: in this case MOOSE will ignore the value you prescribe in controlled_activity_value."}>>> = 'H+'
controlled_activity_value<<<{"description": "Values of the activity or fugacity of the species in controlled_activity_name list. These should always be positive"}>>> = 3.16227E-8 # this is pH=7.5
kinetic_species_name<<<{"description": "Names of the kinetic species given initial values in kinetic_species_initial_value"}>>> = "sulfate_reducer methanogen"
kinetic_species_initial_value<<<{"description": "Initial number of moles, mass or volume (depending on kinetic_species_unit) for each of the species named in kinetic_species_name"}>>> = '1E-6 1E-6'
kinetic_species_unit<<<{"description": "Units of the numerical values given in kinetic_species_initial_value. Moles: mole number. kg: kilograms. g: grams. mg: milligrams. ug: micrograms. cm3: cubic centimeters"}>>> = 'mg mg'
source_species_names<<<{"description": "The name of the species that are added at rates given in source_species_rates. There must be an equal number of these as source_species_rates."}>>> = "H2O Ca++ SO4-- CH3COO- HS- CH4(aq) Fe++"
source_species_rates<<<{"description": "Rates, in mols/time_unit, of addition of the species with names given in source_species_names. A negative value corresponds to removing a species: be careful that you don't cause negative mass problems!"}>>> = "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<<<{"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."}>>> = 1
ramp_max_ionic_strength_subsequent<<<{"description": "The number of iterations over which to progressively increase the maximum ionic strength (from zero to max_ionic_strength) during time-stepping. Unless a great deal occurs in each time step, this parameter can be set quite small"}>>> = 1
execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = ''
solver_info<<<{"description": "Print information (to the console) from the solver including residuals, swaps, etc"}>>> = true
evaluate_kinetic_rates_always<<<{"description": "If true, then, evaluate the kinetic rates at every Newton step during the solve using the current values of molality, activity, etc (ie, implement an implicit solve). If false, then evaluate the kinetic rates using the values of molality, activity, etc, at the start of the current time step (ie, implement an explicit solve)"}>>> = true
adaptive_timestepping<<<{"description": "Use adaptive timestepping at each node in an attempt to ensure convergence of the solver. Setting this parameter to false saves some compute time because copying of datastructures is avoided"}>>> = true
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-14
precision<<<{"description": "Precision for printing values"}>>> = 16
[]
[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"}>>> = 31.536 # 1E-9 mol(acetate)/mg(biomass)/s = 31.536 mol(acetate)/g(biomass)/year
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- SO4--'
promoting_indices<<<{"description": "Indices of the promoting species"}>>> = '1 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 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 200E-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)."}>>> = dissolution
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 # 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<<<{"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
[]
[death_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.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<<<{"description": "Whether the rate should be multiplied by the kinetic_species mass (in grams)"}>>> = true
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)."}>>> = death
eta<<<{"description": "Eta parameter, which appears in |1 - (Q/K)^theta|^eta"}>>> = 0.0
[]
[rate_methanogen]
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"}>>> = "methanogen"
intrinsic_rate_constant<<<{"description": "The intrinsic rate constant for the reaction"}>>> = 63.072 # 2E-9 mol(acetate)/mg(biomass)/s = 63.072 mol(acetate)/g(biomass)/year
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"}>>> = '20E-3'
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
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"}>>> = 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<<<{"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)"}>>> = 24E3
theta<<<{"description": "Theta parameter, which appears in |1 - (Q/K)^theta|^eta"}>>> = 0.5
eta<<<{"description": "Eta parameter, which appears in |1 - (Q/K)^theta|^eta"}>>> = 1
[]
[death_methanogen]
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"}>>> = "methanogen"
intrinsic_rate_constant<<<{"description": "The intrinsic rate constant for the reaction"}>>> = 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<<<{"description": "Whether the rate should be multiplied by the kinetic_species mass (in grams)"}>>> = true
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)."}>>> = death
eta<<<{"description": "Eta parameter, which appears in |1 - (Q/K)^theta|^eta"}>>> = 0.0
[]
[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+ CH3COO- CH4(aq) HS- Ca++ HCO3- SO4-- Fe++"
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 methanogen"
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"}>>> = "*"
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 death_sulfate_reducer rate_methanogen death_methanogen"
[]
[]
[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
type = Transient
[TimeStepper<<<{"href": "../../../syntax/Executioner/TimeStepper/index.html"}>>>]
type = FunctionDT
function = 'min(0.1 * (t + 1), 100)'
[]
end_time = 20000
[]
[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
[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<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
[rate_CH3COO_per_1l_with_source]
type = ParsedAux<<<{"description": "Sets a field variable value to the evaluation of a parsed expression.", "href": "../../../source/auxkernels/ParsedAux.html"}>>>
args = 'rate_CH3COO_per_1l'
variable<<<{"description": "The name of the variable that this object applies to"}>>> = rate_CH3COO_per_1l_with_source
function = 'rate_CH3COO_per_1l + ${rate_CH3COO_diffuse}'
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_begin timestep_end'
[]
[rate_Ca_per_1l_with_source]
type = ParsedAux<<<{"description": "Sets a field variable value to the evaluation of a parsed expression.", "href": "../../../source/auxkernels/ParsedAux.html"}>>>
args = 'rate_Ca_per_1l'
variable<<<{"description": "The name of the variable that this object applies to"}>>> = rate_Ca_per_1l_with_source
function = 'rate_Ca_per_1l + ${rate_Ca_diffuse}'
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_begin timestep_end'
[]
[transported_H2O]
type = GeochemistryQuantityAux<<<{"description": "Extracts information from the Reactor and records it in the AuxVariable", "href": "../../../source/auxkernels/GeochemistryQuantityAux.html"}>>>
variable<<<{"description": "The name of the variable that this object applies to"}>>> = transported_H2O
species<<<{"description": "Species name"}>>> = H2O
quantity<<<{"description": "Type of quantity to output. These are available for non-kinetic species: activity (which equals fugacity for gases); bulk moles (this will be zero if the species is not in the basis); neglog10a = -log10(activity); transported_moles_in_original_basis (will throw an error if species is not in original basis). These are available only for non-kinetic non-minerals: molal = mol(species)/kg(solvent_water); mg_per_kg = mg(species)/kg(solvent_water). These are available only for minerals: free_mg = free mg; free_cm3 = free cubic-centimeters; moles_dumped = moles dumped when special dump and flow-through modes are active. These are available for minerals that host sorbing sites: surface_charge (C/m^2); surface_potential (V). These are available for kinetic species: kinetic_moles; kinetic_additions (-dt * rate = mole increment in kinetic species for this timestep). If quantity=temperature, then the 'species' is ignored and the AuxVariable will record the aqueous solution temperature in degC"}>>> = transported_moles_in_original_basis
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end'
[]
[transported_CH3COO]
type = GeochemistryQuantityAux<<<{"description": "Extracts information from the Reactor and records it in the AuxVariable", "href": "../../../source/auxkernels/GeochemistryQuantityAux.html"}>>>
variable<<<{"description": "The name of the variable that this object applies to"}>>> = transported_CH3COO
species<<<{"description": "Species name"}>>> = "CH3COO-"
quantity<<<{"description": "Type of quantity to output. These are available for non-kinetic species: activity (which equals fugacity for gases); bulk moles (this will be zero if the species is not in the basis); neglog10a = -log10(activity); transported_moles_in_original_basis (will throw an error if species is not in original basis). These are available only for non-kinetic non-minerals: molal = mol(species)/kg(solvent_water); mg_per_kg = mg(species)/kg(solvent_water). These are available only for minerals: free_mg = free mg; free_cm3 = free cubic-centimeters; moles_dumped = moles dumped when special dump and flow-through modes are active. These are available for minerals that host sorbing sites: surface_charge (C/m^2); surface_potential (V). These are available for kinetic species: kinetic_moles; kinetic_additions (-dt * rate = mole increment in kinetic species for this timestep). If quantity=temperature, then the 'species' is ignored and the AuxVariable will record the aqueous solution temperature in degC"}>>> = transported_moles_in_original_basis
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end'
[]
[transported_CH4]
type = GeochemistryQuantityAux<<<{"description": "Extracts information from the Reactor and records it in the AuxVariable", "href": "../../../source/auxkernels/GeochemistryQuantityAux.html"}>>>
variable<<<{"description": "The name of the variable that this object applies to"}>>> = transported_CH4
species<<<{"description": "Species name"}>>> = "CH4(aq)"
quantity<<<{"description": "Type of quantity to output. These are available for non-kinetic species: activity (which equals fugacity for gases); bulk moles (this will be zero if the species is not in the basis); neglog10a = -log10(activity); transported_moles_in_original_basis (will throw an error if species is not in original basis). These are available only for non-kinetic non-minerals: molal = mol(species)/kg(solvent_water); mg_per_kg = mg(species)/kg(solvent_water). These are available only for minerals: free_mg = free mg; free_cm3 = free cubic-centimeters; moles_dumped = moles dumped when special dump and flow-through modes are active. These are available for minerals that host sorbing sites: surface_charge (C/m^2); surface_potential (V). These are available for kinetic species: kinetic_moles; kinetic_additions (-dt * rate = mole increment in kinetic species for this timestep). If quantity=temperature, then the 'species' is ignored and the AuxVariable will record the aqueous solution temperature in degC"}>>> = transported_moles_in_original_basis
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end'
[]
[transported_HS]
type = GeochemistryQuantityAux<<<{"description": "Extracts information from the Reactor and records it in the AuxVariable", "href": "../../../source/auxkernels/GeochemistryQuantityAux.html"}>>>
variable<<<{"description": "The name of the variable that this object applies to"}>>> = transported_HS
species<<<{"description": "Species name"}>>> = "HS-"
quantity<<<{"description": "Type of quantity to output. These are available for non-kinetic species: activity (which equals fugacity for gases); bulk moles (this will be zero if the species is not in the basis); neglog10a = -log10(activity); transported_moles_in_original_basis (will throw an error if species is not in original basis). These are available only for non-kinetic non-minerals: molal = mol(species)/kg(solvent_water); mg_per_kg = mg(species)/kg(solvent_water). These are available only for minerals: free_mg = free mg; free_cm3 = free cubic-centimeters; moles_dumped = moles dumped when special dump and flow-through modes are active. These are available for minerals that host sorbing sites: surface_charge (C/m^2); surface_potential (V). These are available for kinetic species: kinetic_moles; kinetic_additions (-dt * rate = mole increment in kinetic species for this timestep). If quantity=temperature, then the 'species' is ignored and the AuxVariable will record the aqueous solution temperature in degC"}>>> = transported_moles_in_original_basis
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end'
[]
[transported_Ca]
type = GeochemistryQuantityAux<<<{"description": "Extracts information from the Reactor and records it in the AuxVariable", "href": "../../../source/auxkernels/GeochemistryQuantityAux.html"}>>>
variable<<<{"description": "The name of the variable that this object applies to"}>>> = transported_Ca
species<<<{"description": "Species name"}>>> = "Ca++"
quantity<<<{"description": "Type of quantity to output. These are available for non-kinetic species: activity (which equals fugacity for gases); bulk moles (this will be zero if the species is not in the basis); neglog10a = -log10(activity); transported_moles_in_original_basis (will throw an error if species is not in original basis). These are available only for non-kinetic non-minerals: molal = mol(species)/kg(solvent_water); mg_per_kg = mg(species)/kg(solvent_water). These are available only for minerals: free_mg = free mg; free_cm3 = free cubic-centimeters; moles_dumped = moles dumped when special dump and flow-through modes are active. These are available for minerals that host sorbing sites: surface_charge (C/m^2); surface_potential (V). These are available for kinetic species: kinetic_moles; kinetic_additions (-dt * rate = mole increment in kinetic species for this timestep). If quantity=temperature, then the 'species' is ignored and the AuxVariable will record the aqueous solution temperature in degC"}>>> = transported_moles_in_original_basis
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end'
[]
[transported_SO4]
type = GeochemistryQuantityAux<<<{"description": "Extracts information from the Reactor and records it in the AuxVariable", "href": "../../../source/auxkernels/GeochemistryQuantityAux.html"}>>>
variable<<<{"description": "The name of the variable that this object applies to"}>>> = transported_SO4
species<<<{"description": "Species name"}>>> = "SO4--"
quantity<<<{"description": "Type of quantity to output. These are available for non-kinetic species: activity (which equals fugacity for gases); bulk moles (this will be zero if the species is not in the basis); neglog10a = -log10(activity); transported_moles_in_original_basis (will throw an error if species is not in original basis). These are available only for non-kinetic non-minerals: molal = mol(species)/kg(solvent_water); mg_per_kg = mg(species)/kg(solvent_water). These are available only for minerals: free_mg = free mg; free_cm3 = free cubic-centimeters; moles_dumped = moles dumped when special dump and flow-through modes are active. These are available for minerals that host sorbing sites: surface_charge (C/m^2); surface_potential (V). These are available for kinetic species: kinetic_moles; kinetic_additions (-dt * rate = mole increment in kinetic species for this timestep). If quantity=temperature, then the 'species' is ignored and the AuxVariable will record the aqueous solution temperature in degC"}>>> = transported_moles_in_original_basis
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end'
[]
[transported_Fe]
type = GeochemistryQuantityAux<<<{"description": "Extracts information from the Reactor and records it in the AuxVariable", "href": "../../../source/auxkernels/GeochemistryQuantityAux.html"}>>>
variable<<<{"description": "The name of the variable that this object applies to"}>>> = transported_Fe
species<<<{"description": "Species name"}>>> = "Fe++"
quantity<<<{"description": "Type of quantity to output. These are available for non-kinetic species: activity (which equals fugacity for gases); bulk moles (this will be zero if the species is not in the basis); neglog10a = -log10(activity); transported_moles_in_original_basis (will throw an error if species is not in original basis). These are available only for non-kinetic non-minerals: molal = mol(species)/kg(solvent_water); mg_per_kg = mg(species)/kg(solvent_water). These are available only for minerals: free_mg = free mg; free_cm3 = free cubic-centimeters; moles_dumped = moles dumped when special dump and flow-through modes are active. These are available for minerals that host sorbing sites: surface_charge (C/m^2); surface_potential (V). These are available for kinetic species: kinetic_moles; kinetic_additions (-dt * rate = mole increment in kinetic species for this timestep). If quantity=temperature, then the 'species' is ignored and the AuxVariable will record the aqueous solution temperature in degC"}>>> = transported_moles_in_original_basis
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end'
[]
[]
[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
[time]
type = TimePostprocessor<<<{"description": "Reports the current time", "href": "../../../source/postprocessors/TimePostprocessor.html"}>>>
[]
[]
[VectorPostprocessors<<<{"href": "../../../syntax/VectorPostprocessors/index.html"}>>>]
[data]
type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
start_point<<<{"description": "The beginning of the line"}>>> = '0 0 0'
end_point<<<{"description": "The ending of the line"}>>> = '200000 0 0'
num_points<<<{"description": "The number of points to sample along the line"}>>> = 501 # NOTE
sort_by<<<{"description": "What to sort the samples by"}>>> = x
variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = 'transported_CH4 transported_CH3COO transported_SO4 free_mg_sulfate_reducer free_mg_methanogen'
[]
[]
[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
[csv]
type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../source/outputs/CSV.html"}>>>
time_step_interval<<<{"description": "The interval (number of time steps) at which output occurs. Unless explicitly set, the default value of this parameter is set to infinity if the wall_time_interval is explicitly set."}>>> = 10
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = '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<<<{"href": "../../../syntax/Kernels/index.html"}>>>]
[dot_H2O]
type = GeochemistryTimeDerivative<<<{"description": "Kernel describing porosity * d(concentration)/dt, where porosity is an AuxVariable (or just a real number) and concentration is the 'variable' for this Kernel. This Kernel should not be used if porosity is time-dependent. Mass lumping is employed for numerical stability", "href": "../../../source/kernels/GeochemistryTimeDerivative.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_H2O
save_in<<<{"description": "The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)"}>>> = rate_H2O_times_vv
[]
[dot_CH3COO]
type = GeochemistryTimeDerivative<<<{"description": "Kernel describing porosity * d(concentration)/dt, where porosity is an AuxVariable (or just a real number) and concentration is the 'variable' for this Kernel. This Kernel should not be used if porosity is time-dependent. Mass lumping is employed for numerical stability", "href": "../../../source/kernels/GeochemistryTimeDerivative.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_CH3COO
save_in<<<{"description": "The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)"}>>> = rate_CH3COO_times_vv
[]
[dot_CH4]
type = GeochemistryTimeDerivative<<<{"description": "Kernel describing porosity * d(concentration)/dt, where porosity is an AuxVariable (or just a real number) and concentration is the 'variable' for this Kernel. This Kernel should not be used if porosity is time-dependent. Mass lumping is employed for numerical stability", "href": "../../../source/kernels/GeochemistryTimeDerivative.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_CH4
save_in<<<{"description": "The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)"}>>> = rate_CH4_times_vv
[]
[dot_HS]
type = GeochemistryTimeDerivative<<<{"description": "Kernel describing porosity * d(concentration)/dt, where porosity is an AuxVariable (or just a real number) and concentration is the 'variable' for this Kernel. This Kernel should not be used if porosity is time-dependent. Mass lumping is employed for numerical stability", "href": "../../../source/kernels/GeochemistryTimeDerivative.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_HS
save_in<<<{"description": "The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)"}>>> = rate_HS_times_vv
[]
[dot_Ca]
type = GeochemistryTimeDerivative<<<{"description": "Kernel describing porosity * d(concentration)/dt, where porosity is an AuxVariable (or just a real number) and concentration is the 'variable' for this Kernel. This Kernel should not be used if porosity is time-dependent. Mass lumping is employed for numerical stability", "href": "../../../source/kernels/GeochemistryTimeDerivative.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_Ca
save_in<<<{"description": "The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)"}>>> = rate_Ca_times_vv
[]
[dot_SO4]
type = GeochemistryTimeDerivative<<<{"description": "Kernel describing porosity * d(concentration)/dt, where porosity is an AuxVariable (or just a real number) and concentration is the 'variable' for this Kernel. This Kernel should not be used if porosity is time-dependent. Mass lumping is employed for numerical stability", "href": "../../../source/kernels/GeochemistryTimeDerivative.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_SO4
save_in<<<{"description": "The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)"}>>> = rate_SO4_times_vv
[]
[dot_Fe]
type = GeochemistryTimeDerivative<<<{"description": "Kernel describing porosity * d(concentration)/dt, where porosity is an AuxVariable (or just a real number) and concentration is the 'variable' for this Kernel. This Kernel should not be used if porosity is time-dependent. Mass lumping is employed for numerical stability", "href": "../../../source/kernels/GeochemistryTimeDerivative.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_Fe
save_in<<<{"description": "The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)"}>>> = rate_Fe_times_vv
[]
[adv_H2O]
type = ConservativeAdvection<<<{"description": "Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.", "href": "../../../source/kernels/ConservativeAdvection.html"}>>>
velocity = velocity
upwinding_type<<<{"description": "Type of upwinding used. None: Typically results in overshoots and undershoots, but numerical diffusion is minimized. Full: Overshoots and undershoots are avoided, but numerical diffusion is large"}>>> = full
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_H2O
[]
[adv_CH3COO]
type = ConservativeAdvection<<<{"description": "Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.", "href": "../../../source/kernels/ConservativeAdvection.html"}>>>
velocity = velocity
upwinding_type<<<{"description": "Type of upwinding used. None: Typically results in overshoots and undershoots, but numerical diffusion is minimized. Full: Overshoots and undershoots are avoided, but numerical diffusion is large"}>>> = full
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_CH3COO
[]
[adv_CH4]
type = ConservativeAdvection<<<{"description": "Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.", "href": "../../../source/kernels/ConservativeAdvection.html"}>>>
velocity = velocity
upwinding_type<<<{"description": "Type of upwinding used. None: Typically results in overshoots and undershoots, but numerical diffusion is minimized. Full: Overshoots and undershoots are avoided, but numerical diffusion is large"}>>> = full
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_CH4
[]
[adv_HS]
type = ConservativeAdvection<<<{"description": "Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.", "href": "../../../source/kernels/ConservativeAdvection.html"}>>>
velocity = velocity
upwinding_type<<<{"description": "Type of upwinding used. None: Typically results in overshoots and undershoots, but numerical diffusion is minimized. Full: Overshoots and undershoots are avoided, but numerical diffusion is large"}>>> = full
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_HS
[]
[adv_Ca]
type = ConservativeAdvection<<<{"description": "Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.", "href": "../../../source/kernels/ConservativeAdvection.html"}>>>
velocity = velocity
upwinding_type<<<{"description": "Type of upwinding used. None: Typically results in overshoots and undershoots, but numerical diffusion is minimized. Full: Overshoots and undershoots are avoided, but numerical diffusion is large"}>>> = full
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_Ca
[]
[adv_SO4]
type = ConservativeAdvection<<<{"description": "Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.", "href": "../../../source/kernels/ConservativeAdvection.html"}>>>
velocity = velocity
upwinding_type<<<{"description": "Type of upwinding used. None: Typically results in overshoots and undershoots, but numerical diffusion is minimized. Full: Overshoots and undershoots are avoided, but numerical diffusion is large"}>>> = full
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_SO4
[]
[adv_Fe]
type = ConservativeAdvection<<<{"description": "Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.", "href": "../../../source/kernels/ConservativeAdvection.html"}>>>
velocity = velocity
upwinding_type<<<{"description": "Type of upwinding used. None: Typically results in overshoots and undershoots, but numerical diffusion is minimized. Full: Overshoots and undershoots are avoided, but numerical diffusion is large"}>>> = full
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = 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<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
[rate_H2O_auxk]
type = ParsedAux<<<{"description": "Sets a field variable value to the evaluation of a parsed expression.", "href": "../../../source/auxkernels/ParsedAux.html"}>>>
variable<<<{"description": "The name of the variable that this object applies to"}>>> = 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<<<{"href": "../../../syntax/BCs/index.html"}>>>]
[inject_H2O]
type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../source/bcs/DirichletBC.html"}>>>
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'left right'
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = conc_H2O
value<<<{"description": "Value of the BC"}>>> = ${eqm_H2O}
[]
[]At the end of each flow timestep, the geochemistry solver is run:
[MultiApps<<<{"href": "../../../syntax/MultiApps/index.html"}>>>]
[react]
type = TransientMultiApp<<<{"description": "MultiApp for performing coupled simulations with the parent and sub-application both progressing in time.", "href": "../../../source/multiapps/TransientMultiApp.html"}>>>
input_files<<<{"description": "The input file for each App. If this parameter only contains one input file it will be used for all of the Apps. When using 'positions_from_file' it is also admissable to provide one input_file per file."}>>> = bio_zoning_conc.i
clone_parent_mesh<<<{"description": "True to clone parent app mesh and use it for this MultiApp."}>>> = true
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end' # This is critical
[]
[]with the flow rates from the flow solver acting as source terms to the geochemistry simulation:
[Transfers<<<{"href": "../../../syntax/Transfers/index.html"}>>>]
[changes_due_to_flow]
type = MultiAppCopyTransfer<<<{"description": "Copies variables (nonlinear and auxiliary) between multiapps that have identical meshes.", "href": "../../../source/transfers/MultiAppCopyTransfer.html"}>>>
to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = react
source_variable<<<{"description": "The variable to transfer from."}>>> = 'rate_H2O rate_CH3COO rate_CH4 rate_HS rate_Ca rate_SO4 rate_Fe' # change in mole number at every node / dt
variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = '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<<<{"href": "../../../syntax/Transfers/index.html"}>>>]
[transported_moles_from_geochem]
type = MultiAppCopyTransfer<<<{"description": "Copies variables (nonlinear and auxiliary) between multiapps that have identical meshes.", "href": "../../../source/transfers/MultiAppCopyTransfer.html"}>>>
from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = react
source_variable<<<{"description": "The variable to transfer from."}>>> = 'transported_H2O transported_CH3COO transported_CH4 transported_HS transported_Ca transported_SO4 transported_Fe'
variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = '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.
Figure 1: Results of the aquifer-zoning biogeochemistry reactive transport simulation
References
- Craig M. Bethke.
Geochemical and Biogeochemical Reaction Modeling.
Cambridge University Press, 2 edition, 2007.
doi:10.1017/CBO9780511619670.[BibTeX]