The impact of CO fugacity on the solubility of calcite

Section 14.3 of Bethke (2007) and involves studying the solubility of calcite in the presence of a progressively changing CO gas buffer.

Assume that:

  • the bulk composition of Na is 10mmolal;

  • the bulk composition of Cl is 10mmolal;

  • the mineral calcite is used in place of Ca in the basis, and that initially it has a free mole number of 0.01354mol, corresponding to a volume of cm;

  • CO(g) is used in place of H in the basis, and that its fugacity is initially ;

  • charge balance is maintained on HCO.

Then the fugacity of CO(g) is changed from to 1, and the dissolution of calcite, the pH, and species concentrations are all recorded as a function of time.

MOOSE input file

The MOOSE input file includes the GeochemicalModelDefinition which defines the basis species, the equilibrium minerals and equilibrium gases in this case. The piecewise_linear_interpolation flag is set true so that the equilibrium constants and Debye-Huckel parameters are evaluated at the 25C value provided in the database file without any least-squares best-fit smoothing.

[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O H+ Na+ Cl- Ca++ HCO3-"
    equilibrium_minerals = "Calcite"
    equilibrium_gases = "CO2(g)"
    piecewise_linear_interpolation = true # for comparison with GWB
  []
[]
(modules/geochemistry/test/tests/time_dependent_reactions/changing_fugacity_calcite.i)

The TimeDependentReactionSolver specifies

  • the swaps used

  • the initial free molality of Calcite, the initial fugacity of CO(g) and the bulk mole number for the other species

  • the system is closed at by default, meaning that the external agent cannot add or remove Calcite in order to maintain the free mole number of 0.01354 (ie, Calcite can precipitate or dissolve freely)

  • that the activity (fugacity) of CO(g) is controlled by the fug_co2 AuxVariable

[TimeDependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  swap_out_of_basis = "Ca++ H+"
  swap_into_basis = "Calcite CO2(g)"
  charge_balance_species = "HCO3-"
  constraint_species = "H2O              Calcite      CO2(g)        Na+              Cl-              HCO3-"
  constraint_value = "  1.0              0.01354      -3.5          1E-2             1E-2             0"
  constraint_meaning = "kg_solvent_water free_mineral log10fugacity bulk_composition bulk_composition bulk_composition"
  constraint_unit = "   kg               moles        dimensionless moles            moles            moles"
  ramp_max_ionic_strength_initial = 10
  controlled_activity_name = 'CO2(g)'
  controlled_activity_value = fug_co2
  stoichiometric_ionic_str_using_Cl_only = true # for comparison with GWB
  execute_console_output_on = '' # only CSV output required for this example
[]
(modules/geochemistry/test/tests/time_dependent_reactions/changing_fugacity_calcite.i)

The fug_co2 fugacity-controller is

[AuxKernels]
  [fug_co2]
    type = FunctionAux
    variable = fug_co2
    function = '10^(-3.5*(1 - t))'
    execute_on = timestep_begin # so the correct value is provided to the reactor
  []
[]
(modules/geochemistry/test/tests/time_dependent_reactions/changing_fugacity_calcite.i)

with Executioner defining the notion of time:

[Executioner]
  type = Transient
  dt = 0.1
  end_time = 1
[]
(modules/geochemistry/test/tests/time_dependent_reactions/changing_fugacity_calcite.i)

The figures below were produced using a time-step size of 0.01 but the regression test uses a bigger dt for efficiency.

A number of Postprocessors capture information from the AuxVariables automatically included by the TimeDependentReactionSolver:

[Postprocessors]
  [cm3_Calcite]
    type = PointValue
    point = '0 0 0'
    variable = 'free_cm3_Calcite'
  []
  [pH]
    type = PointValue
    point = '0 0 0'
    variable = 'pH'
  []
  [molal_CO2aq]
    type = PointValue
    point = '0 0 0'
    variable = 'molal_CO2(aq)'
  []
  [molal_HCO3-]
    type = PointValue
    point = '0 0 0'
    variable = 'molal_HCO3-'
  []
  [molal_Ca++]
    type = PointValue
    point = '0 0 0'
    variable = 'molal_Ca++'
  []
  [fug_co2]
    type = PointValue
    point = '0 0 0'
    variable = 'activity_CO2(g)'
  []
[]
(modules/geochemistry/test/tests/time_dependent_reactions/changing_fugacity_calcite.i)

GWB input file

The equivalent Geochemists Workbench input file is

# React script that is equivalent to changing_fugacity_calcite.i
Na+ = 1e-2 mol
Cl- = 1e-2 mol
balance on HCO3-
swap Calcite for Ca++
Calcite = 0.01354 free mol
swap CO2(g) for H+
log f CO2(g) = -3.5
epsilon = 1e-13
slide f CO2(g) to 1
go
(modules/geochemistry/test/tests/time_dependent_reactions/changing_fugacity_calcite.rea)

Results

Bethke (2007) presents results in Figures 14.6 and 14.7. The Geochemists Workbench and geochemistry module results are shown below.

Figure 1: Calcite volume as CO fugacity is varied. Compare with Bethke's Figure 14.6

Figure 2: pH of a solution containing calcite as CO fugacity is varied. Compare with Bethke's Figure 14.6

Figure 3: Species concentrations as CO fugacity is varied in a solution containing calcite. Compare with Bethke's Figure 14.7

References

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