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<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [h+]
  []
  [hco3-]
  []
  [ca2+]
  []
  [mg2+]
  []
  [fe2+]
  []
[]
(modules/porous_flow/examples/tutorial/13.i)

The equilibrium reactions are encoded into this Material:

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [equilibrium_massfrac]
    type = PorousFlowMassFractionAqueousEquilibriumChemistry<<<{"description": "This Material forms a std::vector<std::vector ...> of mass-fractions (total concentrations of primary species (m^{3}(primary species)/m^{3}(solution)) and since this is for an aqueous system only, mass-fraction equals volume-fraction) corresponding to an aqueous equilibrium chemistry system.  The first mass fraction is the concentration of the first primary species, etc, and the last mass fraction is the concentration of H2O.", "href": "../../source/materials/PorousFlowMassFractionAqueousEquilibriumChemistry.html"}>>>
    mass_fraction_vars<<<{"description": "List of variables that represent the mass fractions.  For the aqueous phase these are concentrations of the primary species with units m^{3}(chemical)/m^{3}(fluid phase).  For the other phases (if any) these will typically be initialised to zero and will not change throughout the simulation.  Format is 'f_ph0^c0 f_ph0^c1 f_ph0^c2 ... f_ph0^c(N-2) f_ph1^c0 f_ph1^c1 fph1^c2 ... fph1^c(N-2) ... fphP^c0 f_phP^c1 fphP^c2 ... fphP^c(N-2)' where N=number of primary species and P=num_phases, and it is assumed that f_ph^c(N-1)=1-sum(f_ph^c,{c,0,N-2}) so that f_ph^c(N-1) need not be given."}>>> = 'h+ hco3- ca2+ mg2+ fe2+'
    num_reactions<<<{"description": "Number of equations in the system of chemical reactions"}>>> = 5
    equilibrium_constants<<<{"description": "Equilibrium constant for each equation (dimensionless).  If these are temperature dependent AuxVariables, the Jacobian will not be exact"}>>> = 'eqm_k0 eqm_k1 eqm_k2 eqm_k3 eqm_k4'
    primary_activity_coefficients<<<{"description": "Activity coefficients for the primary species (dimensionless) (one for each)"}>>> = '1 1 1 1 1'
    secondary_activity_coefficients<<<{"description": "Activity coefficients for the secondary species (dimensionless) (one for each reaction)"}>>> = '1 1 1 1 1'
    reactions<<<{"description": "A matrix defining the aqueous reactions.  The matrix is entered as a long vector: the first row is entered first, followed by the second row, etc.  There should be num_reactions rows.  All primary species should appear only on the LHS of each reaction (and there should be just one secondary species on the RHS, by definition) so they may have negative coefficients.  Each row should have number of primary_concentrations entries, which are the stoichiometric coefficients.  The first coefficient must always correspond to the first primary species, etc"}>>> = '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:

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [kinetic]
    type = PorousFlowAqueousPreDisChemistry<<<{"description": "This Material forms a std::vector of mineralisation reaction rates (L(precipitate)/L(solution)/s) appropriate to the aqueous precipitation-dissolution system provided.  Note: the PorousFlowTemperature must be measured in Kelvin.", "href": "../../source/materials/PorousFlowAqueousPreDisChemistry.html"}>>>
    primary_concentrations<<<{"description": "List of MOOSE Variables that represent the concentrations of the primary species"}>>> = 'h+ hco3- ca2+ mg2+ fe2+'
    num_reactions<<<{"description": "Number of equations in the system of chemical reactions"}>>> = 1
    equilibrium_constants<<<{"description": "Equilibrium constant for each equation (dimensionless).  If these are temperature dependent AuxVariables, the Jacobian will not be exact"}>>> = kinetic_k
    primary_activity_coefficients<<<{"description": "Activity coefficients for the primary species (dimensionless) (one for each)"}>>> = '1 1 1 1 1'
    reactions<<<{"description": "A matrix defining the aqueous reactions.  The matrix is entered as a long vector: the first row is entered first, followed by the second row, etc.  There should be num_reactions rows.  All primary species should appear only on the LHS of each reaction (and there should be just one secondary species on the RHS, by definition) so they may have negative coefficients.  Each row should have number of primary_concentrations entries, which are the stoichiometric coefficients.  The first coefficient must always correspond to the first primary species, etc"}>>> = '-2 2 1 0.8 0.2'
    specific_reactive_surface_area<<<{"description": "Specific reactive surface area in m^2/(L solution)."}>>> = '1.2E-8'
    kinetic_rate_constant<<<{"description": "Kinetic rate constant in mol/(m^2 s), at the reference temperature (one for each reaction)"}>>> = '3E-4'
    activation_energy<<<{"description": "Activation energy, J/mol (one for each reaction)"}>>> = '1.5e4'
    molar_volume<<<{"description": "Volume occupied by one mole of the secondary species (L(solution)/mol) (one for each reaction)"}>>> = 64365.0
    gas_constant<<<{"description": "Gas constant, in J/(mol K)"}>>> = 8.314
    reference_temperature<<<{"description": "Reference temperature, K"}>>> = 298.15
  []
[]
(modules/porous_flow/examples/tutorial/13.i)

Geometry

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

[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  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:

[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
  [h+_ic]
    type = BoundingBoxIC<<<{"description": "BoundingBoxIC allows setting the initial condition of a value inside and outside of a specified box. The box is aligned with the x, y, z axes", "href": "../../source/ics/BoundingBoxIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = h+
    x1<<<{"description": "The x coordinate of the lower left-hand corner of the box"}>>> = 0.0
    y1<<<{"description": "The y coordinate of the lower left-hand corner of the box"}>>> = 0.0
    x2<<<{"description": "The x coordinate of the upper right-hand corner of the box"}>>> = 1.0e-10
    y2<<<{"description": "The y coordinate of the upper right-hand corner of the box"}>>> = 0.25
    inside<<<{"description": "The value of the variable inside the box"}>>> = 5.0e-2
    outside<<<{"description": "The value of the variable outside the box"}>>> = 1.0e-6
  []
[]
(modules/porous_flow/examples/tutorial/13.i)

and

[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
  [h+_left]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = h+
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 5E-2
  []
[]
(modules/porous_flow/examples/tutorial/13.i)
[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
  [h+_right]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = h+
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    value<<<{"description": "Value of the BC"}>>> = 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<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
  [eqm_k0]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 2.19E6
  []
  [eqm_k1]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 4.73E-11
  []
  [eqm_k2]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0.222
  []
  [eqm_k3]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 1E-2
  []
  [eqm_k4]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 1E-3
  []
  [kinetic_k]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 326.2
  []
  [pressure]
  []
  [dolomite]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [dolomite_initial]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 1E-7
  []
[]
(modules/porous_flow/examples/tutorial/13.i)

and the Material:

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [dolomite_conc]
    type = PorousFlowAqueousPreDisMineral<<<{"description": "This Material forms a std::vector of mineral concentrations (volume-of-mineral/volume-of-material) appropriate to the aqueous precipitation-dissolution system provided.", "href": "../../source/materials/PorousFlowAqueousPreDisMineral.html"}>>>
    initial_concentrations<<<{"description": "Initial concentrations for the mineral species (m^{3}(precipitate)/m^{3}(porous material)).  Default = 0"}>>> = 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.

[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
  [pressure_ic]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pressure
    function<<<{"description": "The initial condition 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

Two of the primary species concentrations at the end of the simulation.

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

The precipitated dolomite concentration at the end of the simulation.

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

Start | Previous |