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]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../database/selenate_sorption.json"
    basis_species = "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]
  model_definition = definition
  charge_balance_species = "Na+"
  constraint_species = "H2O              H+            Na+              SorbingSite      SeO4--"
  constraint_value = "  1.0              -7.5          10E-6            310E-9           5E-6"
  constraint_meaning = "kg_solvent_water log10activity bulk_composition bulk_composition free_concentration"
  constraint_unit = "   kg               dimensionless moles            moles            molal"
  ramp_max_ionic_strength_initial = 0 # not needed in this simple problem
  execute_console_output_on = '' # 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]
  [mol_sorbed_selenate_per_g_dry_soil]
  []
[]
[AuxKernels]
  [mol_sorbed_selenate_per_g_dry_soil]
    type = ParsedAux
    coupled_variables = molal_SorbedSelenate
    expression = 'molal_SorbedSelenate / 500.0'
    variable = mol_sorbed_selenate_per_g_dry_soil
  []
[]
[Postprocessors]
  [mol_sorbed_selenate_per_g_dry_soil]
    type = PointValue
    point = '0 0 0'
    variable = mol_sorbed_selenate_per_g_dry_soil
  []
[]

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