Porous Flow Tutorial Page 07. A chemically-reactive tracer, and porosity and permeability changes
In this Page, the model from Page 06 is enhanced so that its tracer is chemically reactive. PorousFlow has the ability to simulate aqueous equilibrium chemistry and aqueous precipitation-dissolution chemistry.
As an illustration, in this Page, the tracer is assumed to precipitate according to
(1)
where is the tracer concentration (m(tracer)/m(solution)) is a mineral (m(mineral)/m(porous-material)). This precipitation causes porosity changes and permeability changes.
The input file of Page 06 is modified to include chemistry. The first modification is to the PorousFlowFullySaturated so that it understands there is one chemical reaction:
[PorousFlowFullySaturated<<<{"href": "../../syntax/PorousFlowFullySaturated/index.html"}>>>]
porepressure<<<{"description": "The name of the porepressure variable"}>>> = porepressure
coupling_type<<<{"description": "The type of simulation. For simulations involving Mechanical deformations, you will need to supply the correct Biot coefficient. For simulations involving Thermal flows, you will need an associated ConstantThermalExpansionCoefficient Material"}>>> = Hydro
gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = the_simple_fluid
mass_fraction_vars<<<{"description": "List of variables that represent the mass fractions. With only one fluid component, this may be left empty. With N fluid components, the format is 'f_0 f_1 f_2 ... f_(N-1)'. That is, the N^th component need not be specified because f_N = 1 - (f_0 + f_1 + ... + f_(N-1)). It is best numerically to choose the N-1 mass fraction variables so that they represent the fluid components with small concentrations. This Action will associated the i^th mass fraction variable to the equation for the i^th fluid component, and the pressure variable to the N^th fluid component."}>>> = tracer_concentration
number_aqueous_kinetic<<<{"description": "The number of secondary species in the aqueous-kinetic reaction system involved in precipitation and dissolution. (Leave as zero if the simulation does not involve chemistry)"}>>> = 1
temperature<<<{"description": "For isothermal simulations, this is the temperature at which fluid properties (and stress-free strains) are evaluated at. Otherwise, this is the name of the temperature variable. Units = Kelvin"}>>> = 283.0
stabilization<<<{"description": "Numerical stabilization used. 'Full' means full upwinding. 'KT' means FEM-TVD stabilization of Kuzmin-Turek"}>>> = none # Note to reader: try this with other stabilization and compare the results
[]Setting the temperature here allows control of the chemical reaction rate (which is temperature dependent: see chemical reactions).
As precipitation or dissolution occurs via Eq. (1), the tracer concentration gets decreased or increased. Thus, a new Kernel must be added, which is a PorousFlowPreDis Kernel:
[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
[precipitation_dissolution]
type = PorousFlowPreDis<<<{"description": "Precipitation-dissolution of chemical species", "href": "../../source/kernels/PorousFlowPreDis.html"}>>>
mineral_density<<<{"description": "Density (kg(precipitate)/m^3(precipitate)) of each secondary species in the aqueous precipitation-dissolution reaction system"}>>> = 1000.0
stoichiometry<<<{"description": "A vector of stoichiometric coefficients for the primary species that is the Variable of this Kernel: one for each precipitation-dissolution reaction (these are one columns of the 'reactions' matrix)"}>>> = 1
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tracer_concentration
[]
[]The chemical reaction rate is computed by a PorousFlowAqueousPreDisChemistry Material:
[precipitation_dissolution_mat]
type = PorousFlowAqueousPreDisChemistry
reference_temperature = 283.0
activation_energy = 1 # irrelevant because T=Tref
equilibrium_constants = eqm_k # equilibrium tracer concentration
kinetic_rate_constant = 1E-8
molar_volume = 1
num_reactions = 1
primary_activity_coefficients = 1
primary_concentrations = tracer_concentration
reactions = 1
specific_reactive_surface_area = 1
[]All the parameters are fully explained in the chemical reactions module. Briefly, in this case, because the reference_temperature equals the temperature specified in the PorousFlowFullySaturated, there is no temperature dependence of the reaction rate, so it is just the product of the kinetic_rate_constant (mol.m.s), the specific_reactive_surface_area (1m.L), the molar volume (1L.mol) and , where is the equilibrium constant (0.1): This Material feeds its rate into the PorousFlowPreDis Kernel (to alter the tracer concentration), as well as into a PorousFlowAqueousPreDisMineral Material:
[mineral_concentration]
type = PorousFlowAqueousPreDisMineral
The reason for this Material is that we can build an AuxVariable (below) to record the concentration of the precipitated mineral, but also because the mineral concentration enters into the porosity calculation. The porosity material is
[porosity_mat]
type = PorousFlowPorosity
porosity_zero = 0.1
chemical = true
initial_mineral_concentrations = initial_and_reference_conc
reference_chemistry = initial_and_reference_conc
[]In PorousFlow, the permeability may depend on the porosity via the KozenyCarman relationship. In this case, the Materials look like:
[permeability_aquifer]
type = PorousFlowPermeabilityKozenyCarman
block = aquifer
k0 = 1E-14
m = 2
n = 3
phi0 = 0.1
poroperm_function = kozeny_carman_phi0
[]
[permeability_caps]
type = PorousFlowPermeabilityKozenyCarman
block = caps
k0 = 1E-15
k_anisotropy = '1 0 0 0 1 0 0 0 0.1'
m = 2
n = 3
phi0 = 0.1
poroperm_function = kozeny_carman_phi0
[]Instead of using a set of preset DirichletBC to control the physics at the injection area, tracer is simply injected using a PorousFlowSink (see also boundary conditions for a detailed discussion). A fixed rate of kg.s.m is used:
[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
[constant_injection_of_tracer]
type = PorousFlowSink<<<{"description": "Applies a flux sink to a boundary.", "href": "../../source/bcs/PorousFlowSink.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tracer_concentration
flux_function<<<{"description": "The flux. The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE. The functional form is useful for spatially or temporally varying sinks. Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = -5E-3
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
[]
[constant_outer_porepressure]
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"}>>> = porepressure
value<<<{"description": "Value of the BC"}>>> = 0
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = rmax
[]
[]
[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
[the_simple_fluid]
type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../source/fluidproperties/SimpleFluidProperties.html"}>>>
bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2E9
viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1.0E-3
density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000.0
[]
[]
[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
[porosity_mat]
type = PorousFlowPorosity<<<{"description": "This Material calculates the porosity PorousFlow simulations", "href": "../../source/materials/PorousFlowPorosity.html"}>>>
porosity_zero<<<{"description": "The porosity at zero volumetric strain and reference temperature and reference effective porepressure and reference chemistry. This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)"}>>> = 0.1
chemical<<<{"description": "If true, porosity will be a function of precipitate"}>>> = true
initial_mineral_concentrations<<<{"description": "Initial mineral concentrations (m^3(precipitate)/m^3(porous material)), entered as a vector (one value per mineral). (Only used if chemical=true)"}>>> = initial_and_reference_conc
reference_chemistry<<<{"description": "Reference values of the solid mineral concentrations (m^3(precipitate)/m^3(porous material)), entered as a vector (one value per mineral). (Only used if chemical=true)"}>>> = initial_and_reference_conc
[]
[permeability_aquifer]
type = PorousFlowPermeabilityKozenyCarman<<<{"description": "This Material calculates the permeability tensor from a form of the Kozeny-Carman equation based on the spatially constant initial permeability and porosity or grain size.", "href": "../../source/materials/PorousFlowPermeabilityKozenyCarman.html"}>>>
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = aquifer
k0<<<{"description": "The permeability scalar value (usually in m^2) at the reference porosity, required for kozeny_carman_phi0"}>>> = 1E-14
m<<<{"description": "(1-porosity) exponent"}>>> = 2
n<<<{"description": "Porosity exponent"}>>> = 3
phi0<<<{"description": "The reference porosity, required for kozeny_carman_phi0"}>>> = 0.1
poroperm_function<<<{"description": "Function relating porosity and permeability. The options are: kozeny_carman_fd2 = f d^2 phi^n/(1-phi)^m (where phi is porosity, f is a scalar constant with typical values 0.01-0.001, and d is grain size). kozeny_carman_phi0 = k0 (1-phi0)^m/phi0^n * phi^n/(1-phi)^m (where phi is porosity, and k0 is the permeability at porosity phi0) kozeny_carman_A = A for directly supplying the permeability multiplying factor."}>>> = kozeny_carman_phi0
[]
[permeability_caps]
type = PorousFlowPermeabilityKozenyCarman<<<{"description": "This Material calculates the permeability tensor from a form of the Kozeny-Carman equation based on the spatially constant initial permeability and porosity or grain size.", "href": "../../source/materials/PorousFlowPermeabilityKozenyCarman.html"}>>>
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = caps
k0<<<{"description": "The permeability scalar value (usually in m^2) at the reference porosity, required for kozeny_carman_phi0"}>>> = 1E-15
k_anisotropy<<<{"description": "A tensor to multiply the calculated scalar permeability, in order to obtain anisotropy if required. Defaults to isotropic permeability if not specified."}>>> = '1 0 0 0 1 0 0 0 0.1'
m<<<{"description": "(1-porosity) exponent"}>>> = 2
n<<<{"description": "Porosity exponent"}>>> = 3
phi0<<<{"description": "The reference porosity, required for kozeny_carman_phi0"}>>> = 0.1
poroperm_function<<<{"description": "Function relating porosity and permeability. The options are: kozeny_carman_fd2 = f d^2 phi^n/(1-phi)^m (where phi is porosity, f is a scalar constant with typical values 0.01-0.001, and d is grain size). kozeny_carman_phi0 = k0 (1-phi0)^m/phi0^n * phi^n/(1-phi)^m (where phi is porosity, and k0 is the permeability at porosity phi0) kozeny_carman_A = A for directly supplying the permeability multiplying factor."}>>> = kozeny_carman_phi0
[]
[precipitation_dissolution_mat]
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"}>>>
reference_temperature<<<{"description": "Reference temperature, K"}>>> = 283.0
activation_energy<<<{"description": "Activation energy, J/mol (one for each reaction)"}>>> = 1 # irrelevant because T=Tref
equilibrium_constants<<<{"description": "Equilibrium constant for each equation (dimensionless). If these are temperature dependent AuxVariables, the Jacobian will not be exact"}>>> = eqm_k # equilibrium tracer concentration
kinetic_rate_constant<<<{"description": "Kinetic rate constant in mol/(m^2 s), at the reference temperature (one for each reaction)"}>>> = 1E-8
molar_volume<<<{"description": "Volume occupied by one mole of the secondary species (L(solution)/mol) (one for each reaction)"}>>> = 1
num_reactions<<<{"description": "Number of equations in the system of chemical reactions"}>>> = 1
primary_activity_coefficients<<<{"description": "Activity coefficients for the primary species (dimensionless) (one for each)"}>>> = 1
primary_concentrations<<<{"description": "List of MOOSE Variables that represent the concentrations of the primary species"}>>> = tracer_concentration
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
specific_reactive_surface_area<<<{"description": "Specific reactive surface area in m^2/(L solution)."}>>> = 1
[]
[mineral_concentration]
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"}>>>
[]
[]
[Preconditioning<<<{"href": "../../syntax/Preconditioning/index.html"}>>>]
active<<<{"description": "If specified only the blocks named will be visited and made active"}>>> = basic
[basic]
type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = ' lu mumps'
[]
[]
[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-10
[]
[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
[]Finally, a set of AuxVariables (some employing the handy PorousFlowPropertyAux) is used to record useful information:
[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
[eqm_k]
initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0.1
[]
[mineral_conc]
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
[]
[initial_and_reference_conc]
initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0
[]
[porosity]
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
[]
[permeability]
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
[]
[]
[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
[mineral_conc]
type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = mineral_concentration
mineral_species<<<{"description": "The mineral chemical species number"}>>> = 0
variable<<<{"description": "The name of the variable that this object applies to"}>>> = mineral_conc
[]
[porosity]
type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = porosity
variable<<<{"description": "The name of the variable that this object applies to"}>>> = porosity
[]
[permeability]
type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = permeability
column<<<{"description": "Column of permeability tensor to output"}>>> = 0
row<<<{"description": "Row of permeability tensor to output"}>>> = 0
variable<<<{"description": "The name of the variable that this object applies to"}>>> = permeability
[]
[]An animation of the results is shown in Figure 1.
Figure 1: Tracer concentration, porepressure, mineral concentration and permeability evolution in the borehole-aquifer-caprock system.
The simulation may be promoted to a full THMC simulation using the approach used in Page 03 and Page 04. The arguments made about scaling the variables must be modified to take into account the fluid density appearing in the fluid equation (see Page 06) so the scaling will be approximately for the temperature and for the displacement variables. The porosity may depend on porepressure, temperature, volumetric strain and chemistry.