Sorption of selenate in loamy soil

This example closely follows Section 9.6 of Bethke (2007).

We explore the sorption of selenate (SeO) as predicted by the Langmuir approach. It is assumed

• the sample consists of 1kg of water and 500g of dry soil (corresponding to 189cm of dry soil);

• the number of moles of sorbing sites per gram of dry soil is mol.g. This sets the total mole number of sorbing sites, mol which appears in the MOOSE input file;

• the equilibrium constant for Langmuir sorption is ;

• the pH is 7.5;

• charge balance is enforced on Na.

Selenate is a redox couple in the database. Decoupling it, the basis is (H20, H+, Na+, selenate, SorbingSite). The sorbing reaction, involving selenate, SorbingSite and sorbedSelenate, is most conveniently introduced into MOOSE by an additional database file that includes the aforementioned equilibrium constant. The GeochemicalModelDefinition therefore reads:

[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/selenate_sorption.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+ Na+ SeO4-- SorbingSite"
  []
[]

The TimeIndependentReactionSolver sets:

• the extent of the system to 1kg of solvent water;

• the pH (via the H activity);

• a rough initial condition for the charge-balance species Na;

• the total mole number of sorbing sites to mol;

• the free molality of selenate, which is the independent variable in this problem.

[TimeIndependentReactionSolver<<<{"href": "../../../syntax/TimeIndependentReactionSolver/index.html"}>>>]
  model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object (you must create this UserObject yourself)"}>>> = definition
  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."}>>> = "Na+"
  constraint_species<<<{"description": "Names of the species that have their values fixed to constraint_value with meaning constraint_meaning.  All basis species (after swap_into_basis and swap_out_of_basis) must be provided with exactly one constraint.  These constraints are used to compute the configuration during the initial problem setup, and in time-dependent simulations they may be modified as time progresses."}>>> = "H2O              H+            Na+              SorbingSite      SeO4--"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              -7.5          10E-6            310E-9           5E-6"
  constraint_meaning<<<{"description": "Meanings of the numerical values given in constraint_value.  kg_solvent_water: can only be applied to H2O and units must be kg.  bulk_composition: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species but not in kinetic species, and units must be moles or mass (kg, g, etc).  bulk_composition_with_kinetic: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species and in kinetic species, and units must be moles or mass (kg, g, etc).  free_concentration: can be applied to all basis species that are not gas and not H2O and not mineral, and represents the total amount of the basis species existing freely (not as secondary species) within the solution, and units must be molal or mass_per_kg_solvent.  free_mineral: can be applied to all mineral basis species, and represents the total amount of the mineral existing freely (precipitated) within the solution, and units must be moles, mass or cm3.  activity and log10activity: can be applied to basis species that are not gas and not mineral and not sorbing sites, and represents the activity of the basis species (recall pH = -log10activity), and units must be dimensionless.  fugacity and log10fugacity: can be applied to gases, and units must be dimensionless"}>>> = "kg_solvent_water log10activity bulk_composition bulk_composition free_concentration"
  constraint_unit<<<{"description": "Units of the numerical values given in constraint_value.  Dimensionless: should only be used for activity or fugacity constraints.  Moles: mole number.  Molal: moles per kg solvent water.  kg: kilograms.  g: grams.  mg: milligrams.  ug: micrograms.  kg_per_kg_solvent: kilograms per kg solvent water.  g_per_kg_solvent: grams per kg solvent water.  mg_per_kg_solvent: milligrams per kg solvent water.  ug_per_kg_solvent: micrograms per kg solvent water.  cm3: cubic centimeters"}>>> = "   kg               dimensionless moles            moles            molal"
  ramp_max_ionic_strength_initial<<<{"description": "The number of iterations over which to progressively increase the maximum ionic strength (from zero to max_ionic_strength) during the initial equilibration.  Increasing this can help in convergence of the Newton process, at the cost of spending more time finding the aqueous configuration."}>>> = 0 # not needed in this simple problem
  execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = '' # only CSV output for this problem
[]

The Langmuir approach uses a reaction of the form with an equilibrium constant of in this case. Assuming all activity coefficients are equal to 1.0 means the molalities, , are related via (1) The equation for the bulk composition (mole number) of SorbingSite is where is the mass of the solvent water. Substituting this into the Eq. (1) yields where .

The analysis of the preceding paragraph fairly accurately describes the current situation because the ionic strength is close to zero, so the activities are all close to unity. The MOOSE input file uses the molal_SorbedSelenate AuxVariable produced by the TimeIndependentReactionSolver to obtain the moles of sorbed selenate per gram of dry soil:

[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
  [mol_sorbed_selenate_per_g_dry_soil]
  []
[]
[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [mol_sorbed_selenate_per_g_dry_soil]
    type = ParsedAux<<<{"description": "Sets a field variable value to the evaluation of a parsed expression.", "href": "../../../source/auxkernels/ParsedAux.html"}>>>
    coupled_variables<<<{"description": "Vector of coupled variable names"}>>> = molal_SorbedSelenate
    expression<<<{"description": "Parsed function expression to compute"}>>> = 'molal_SorbedSelenate / 500.0'
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = mol_sorbed_selenate_per_g_dry_soil
  []
[]
[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [mol_sorbed_selenate_per_g_dry_soil]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = mol_sorbed_selenate_per_g_dry_soil
  []
[]

Running the simulation with different free_molality for the selenate produces the results shown in Table 1.