Start | Previous |

Porous Flow Tutorial Page 13. More elaborate chemistry

A very simple chemical system was built in Page 07. The reader is encouraged to consult that page before moving to the more elaborate situation described here. This page:

  • does not use an Action to describe the Kernels, Materials, etc. This makes the input file quite long, but perhaps more easy to extend to multi-phase situations, different boundary conditions, etc. An Action could easily be used instead, by copying the chemistry from this tutorial to that on Page 07.

  • builds a fully-saturated aqueous chemical system that could be used to describe dolomite precipitation and dissolution.

This page illustrates that it is possible to build quite complicated chemical systems within PorousFlow. However, geochemists will recognize these are still unrealistically simple, since they contain only a few species, the activity coefficients are all unity, the equilibrium constants are not easily temperature-dependent, etc. If you require state-of-the-art geochemical modelling capability, please use MOOSE's Geochemistry module.

The equilibrium system

The equilibrium system has:

  • 5 primary species, which are H, HCO, Ca, Mg, Fe.

  • 5 secondary species, which are CO(aq), CO, CaHCO, MgHCO, FeHCO.

The equations are Some of these equilibrium constants have been chosen rather arbitrarily.

The primary species are represented as PorousFlow variables:

[Variables]
  [h+]
  []
  [hco3-]
  []
  [ca2+]
  []
  [mg2+]
  []
  [fe2+]
  []
[]
(modules/porous_flow/examples/tutorial/13.i)

The equilibrium reactions are encoded into this Material:

  [equilibrium_massfrac]
    type = PorousFlowMassFractionAqueousEquilibriumChemistry
    mass_fraction_vars = 'h+ hco3- ca2+ mg2+ fe2+'
    num_reactions = 5
    equilibrium_constants = 'eqm_k0 eqm_k1 eqm_k2 eqm_k3 eqm_k4'
    primary_activity_coefficients = '1 1 1 1 1'
    secondary_activity_coefficients = '1 1 1 1 1'
    reactions = '1 1 0 0 0
                -1 1 0 0 0
                 0 1 1 0 0
                 0 1 0 1 0
                 0 1 0 0 1'
  []
(modules/porous_flow/examples/tutorial/13.i)

The kinetic system

This is with the following parameters:

  • molar volume 64365.0L(solution)/mol,

  • mineral density 2875.0kg(precipitate)/m(precipitate)

  • equilibrium constant ,

  • specific reactive surface area m/L,

  • kinetic rate constant mol/m/s,

  • activation energy J/mol,

  • the primary activity coefficients are all unity,

  • and the and exponents are also unity.

Some of these quantities have been chosen rather arbitrarily. This kinetic system is encoded in the input file as:

  [kinetic]
    type = PorousFlowAqueousPreDisChemistry
    primary_concentrations = 'h+ hco3- ca2+ mg2+ fe2+'
    num_reactions = 1
    equilibrium_constants = kinetic_k
    primary_activity_coefficients = '1 1 1 1 1'
    reactions = '-2 2 1 0.8 0.2'
    specific_reactive_surface_area = '1.2E-8'
    kinetic_rate_constant = '3E-4'
    activation_energy = '1.5e4'
    molar_volume = 64365.0
    gas_constant = 8.314
    reference_temperature = 298.15
  []
(modules/porous_flow/examples/tutorial/13.i)

Geometry

The model is just a 1D line, extending between and .

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 100
  xmax = 1
[]
(modules/porous_flow/examples/tutorial/13.i)

The initial and boundary conditions

Primary variables

Each of the primary variables are initialised to have concentration m(species)/m(solution) everywhere in the domain except for at the left-hand side () where they have concentration . The boundary conditions are to fix these values at the left and right sides of the domain. For instance:

  [h+_ic]
    type = BoundingBoxIC
    variable = h+
    x1 = 0.0
    y1 = 0.0
    x2 = 1.0e-10
    y2 = 0.25
    inside = 5.0e-2
    outside = 1.0e-6
  []
(modules/porous_flow/examples/tutorial/13.i)

and

  [h+_left]
    type = DirichletBC
    variable = h+
    boundary = left
    value = 5E-2
  []
(modules/porous_flow/examples/tutorial/13.i)
  [h+_right]
    type = DirichletBC
    variable = h+
    boundary = right
    value = 1e-6
  []
(modules/porous_flow/examples/tutorial/13.i)

Please remember that boundary conditions in PorousFlow are usually more complicated than setting Dirichlet or Preset boundary conditions: see boundary conditions. Looking at the results below you can clearly see the effect of the naive boundary conditions placed on the right-hand side.

Dolomite

The initial condition for dolomite is m(precipitate)/m(porous material). This is implemented in the Auxiliary system by

[AuxVariables]
  [eqm_k0]
    initial_condition = 2.19E6
  []
  [eqm_k1]
    initial_condition = 4.73E-11
  []
  [eqm_k2]
    initial_condition = 0.222
  []
  [eqm_k3]
    initial_condition = 1E-2
  []
  [eqm_k4]
    initial_condition = 1E-3
  []
  [kinetic_k]
    initial_condition = 326.2
  []
  [pressure]
  []
  [dolomite]
    family = MONOMIAL
    order = CONSTANT
  []
  [dolomite_initial]
    initial_condition = 1E-7
  []
[]

[AuxKernels]
  [dolomite]
    type = PorousFlowPropertyAux
    property = mineral_concentration
    mineral_species = 0
    variable = dolomite
  []
[]
(modules/porous_flow/examples/tutorial/13.i)

and the Material:

  [dolomite_conc]
    type = PorousFlowAqueousPreDisMineral
    initial_concentrations = dolomite_initial
  []
(modules/porous_flow/examples/tutorial/13.i)

Given the above equilibrium constant and concentrations of the primary species, the dolomite immediately begins to dissolve into solution.

Porepressure

The porepressure is fixed to have gradient Pa/m.

  [pressure_ic]
    type = FunctionIC
    variable = pressure
    function = '(1 - x) * 1E6'
  []
(modules/porous_flow/examples/tutorial/13.i)

With a permeability of m and a fluid viscosity of Pa.s, the Darcy velocity is m/s. The porosity is held fixed at 0.2.

Results

Figure 1: Two of the primary species concentrations at the end of the simulation.

Figure 2: The precipitated dolomite concentration at the end of the simulation.

Start | Previous |