Exploring the impact of pH on sorption

Section 14.3 of Bethke (2007) describes how to explore the impact of pH on sorption.

This page builds upon the surface complexation example in which sorption to the ferric hydroxide mineral Fe(OH) was studied. In that example, agreement between Bethke (2007), the Geochemists Workbench software, and the geochemistry module was demonstrated at pH 4 and 8. On this page, the pH is varied between 4 and 12, and amount of sorption of each species is plotted.

The surface complexation example included lead, mercury and iron complexes, but for simplicity this example only include iron complexes. The model definition therefore reads:

[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../database/ferric_hydroxide_sorption.json"
    basis_species = "H2O H+ Na+ Cl- Fe+++ >(s)FeOH >(w)FeOH"
    equilibrium_minerals = "Fe(OH)3(ppd)"
    piecewise_linear_interpolation = true # for comparison with GWB
  []
[]
(modules/geochemistry/test/tests/time_dependent_reactions/changing_pH_ferric_hydroxide.i)

The TimeDependentReactionSolver defines the initial composition as well as how to vary the pH with time via the controlled_activity inputs:

[TimeDependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  swap_out_of_basis = "Fe+++"
  swap_into_basis = "Fe(OH)3(ppd)"
  charge_balance_species = "Cl-"
  constraint_species = "H2O              H+            Na+              Cl-              Fe(OH)3(ppd) >(s)FeOH         >(w)FeOH"
  constraint_value = "  1.0              -4            0.1              0.1              9.3573E-3    4.6786E-5        1.87145E-3"
  constraint_meaning = "kg_solvent_water log10activity bulk_composition bulk_composition free_mineral bulk_composition bulk_composition"
  constraint_unit = "kg               dimensionless moles            moles            moles        moles            moles"
  controlled_activity_name = "H+"
  controlled_activity_value = set_aH
  ramp_max_ionic_strength_initial = 0 # not needed in this simple problem
  stoichiometric_ionic_str_using_Cl_only = true # for comparison with GWB
  abs_tol = 1E-14
  execute_console_output_on = '' # only CSV output needed for this example
[]
(modules/geochemistry/test/tests/time_dependent_reactions/changing_pH_ferric_hydroxide.i)

The set_aH value is defined via a FunctionAux in the following way:

[AuxVariables]
  [set_aH]
  []
[]
[AuxKernels]
  [set_aH]
    type = FunctionAux
    variable = set_aH
    function = '10^(-4-t)'
    execute_on = timestep_begin # so the correct value is provided to the reactor
  []
[]
(modules/geochemistry/test/tests/time_dependent_reactions/changing_pH_ferric_hydroxide.i)

and the timestepping is defined in the executioner:

[Executioner]
  type = Transient
  start_time = -0.25
  dt = 0.25
  end_time = 8
[]
(modules/geochemistry/test/tests/time_dependent_reactions/changing_pH_ferric_hydroxide.i)

In this situation, the start_time is before when the system closes and the time-dependency begins, just so the results at can be simply recorded. Finally, the quantities of interest are recorded into Postprocessors using the AuxVariables that are automatically included in the simulation by the TimeDependentReactionSolver

[Postprocessors]
  [pH]
    type = PointValue
    point = '0 0 0'
    variable = 'pH'
  []
  [molal_>wFeOH2+]
    type = PointValue
    point = '0 0 0'
    variable = 'molal_>(w)FeOH2+'
  []
  [molal_>wFeOH]
    type = PointValue
    point = '0 0 0'
    variable = 'molal_>(w)FeOH'
  []
  [molal_>wFeO-]
    type = PointValue
    point = '0 0 0'
    variable = 'molal_>(w)FeO-'
  []
  [molal_>sFeOH2+]
    type = PointValue
    point = '0 0 0'
    variable = 'molal_>(s)FeOH2+'
  []
  [molal_>sFeOH]
    type = PointValue
    point = '0 0 0'
    variable = 'molal_>(s)FeOH'
  []
  [molal_>sFeO-]
    type = PointValue
    point = '0 0 0'
    variable = 'molal_>(s)FeO-'
  []
  [potential]
    type = PointValue
    point = '0 0 0'
    variable = 'surface_potential_Fe(OH)3(ppd)'
  []
[]
(modules/geochemistry/test/tests/time_dependent_reactions/changing_pH_ferric_hydroxide.i)

Bethke (2007) presents results in Figures 14.8 and 14.9 (bold line only, since the fine line shows a different type of aqueous solution). The results are faithfully reproduced by the geochemistry module as shown in the figures below.

Figure 1: Concentrations of sites on a ferric oxide surface. Compare with Bethke's Figure 14.8

Figure 2: Surface potential of ferric oxide. Compare with Bethke's Figure 14.9

References

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