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
(modules/geochemistry/test/tests/time_dependent_reactions/changing_pH_ferric_hydroxide.i)
# Sorption onto FerricHydroxide along with changing pH
[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
[]
[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
[]
[]
[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)'
[]
[]
[Executioner]
type = Transient
start_time = -0.25
dt = 0.25
end_time = 8
[]
[Outputs]
csv = true
[]
[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)
# Sorption onto FerricHydroxide along with changing pH
[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
[]
[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
[]
[]
[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)'
[]
[]
[Executioner]
type = Transient
start_time = -0.25
dt = 0.25
end_time = 8
[]
[Outputs]
csv = true
[]
[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)
# Sorption onto FerricHydroxide along with changing pH
[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
[]
[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
[]
[]
[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)'
[]
[]
[Executioner]
type = Transient
start_time = -0.25
dt = 0.25
end_time = 8
[]
[Outputs]
csv = true
[]
[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)
# Sorption onto FerricHydroxide along with changing pH
[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
[]
[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
[]
[]
[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)'
[]
[]
[Executioner]
type = Transient
start_time = -0.25
dt = 0.25
end_time = 8
[]
[Outputs]
csv = true
[]
[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)
# Sorption onto FerricHydroxide along with changing pH
[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
[]
[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
[]
[]
[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)'
[]
[]
[Executioner]
type = Transient
start_time = -0.25
dt = 0.25
end_time = 8
[]
[Outputs]
csv = true
[]
[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
[]
[]