Start | Previous | Next

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]
  porepressure = porepressure
  coupling_type = Hydro
  gravity = '0 0 0'
  fp = the_simple_fluid
  mass_fraction_vars = tracer_concentration
  number_aqueous_kinetic = 1
  temperature = 283.0
[]
(modules/porous_flow/examples/tutorial/07.i)

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]
  [./precipitation_dissolution]
    type = PorousFlowPreDis
    mineral_density = 1000.0
    stoichiometry = 1
    variable = tracer_concentration
  [../]
[]
(modules/porous_flow/examples/tutorial/07.i)

The chemical reaction rate is computed by a PorousFlowAqueousPreDisChemistry Material:

# Darcy flow with a tracer that precipitates causing mineralisation and porosity changes and permeability changes
[Mesh]
  type = AnnularMesh
  dim = 2
  nr = 10
  rmin = 1.0
  rmax = 10
  growth_r = 1.4
  nt = 4
  tmin = 0
  tmax = 1.57079632679
[]

[MeshModifiers]
  [./make3D]
    type = MeshExtruder
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
  [../]
  [./shift_down]
    type = Transform
    transform = TRANSLATE
    vector_value = '0 0 -6'
    depends_on = make3D
  [../]
  [./aquifer]
    type = SubdomainBoundingBox
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    depends_on = shift_down
  [../]
  [./injection_area]
    type = ParsedAddSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    depends_on = 'aquifer'
  [../]
  [./rename]
    type = RenameBlock
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    depends_on = 'injection_area'
  [../]
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [./porepressure]
  [../]
  [./tracer_concentration]
  [../]
[]

[PorousFlowFullySaturated]
  porepressure = porepressure
  coupling_type = Hydro
  gravity = '0 0 0'
  fp = the_simple_fluid
  mass_fraction_vars = tracer_concentration
  number_aqueous_kinetic = 1
  temperature = 283.0
[]

[AuxVariables]
  [./eqm_k]
    initial_condition = 0.1
  [../]
  [./mineral_conc]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./initial_and_reference_conc]
    initial_condition = 0
  [../]
  [./porosity]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./permeability]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[AuxKernels]
  [./mineral_conc]
    type = PorousFlowPropertyAux
    property = mineral_concentration
    mineral_species = 0
    variable = mineral_conc
  [../]
  [./porosity]
    type = PorousFlowPropertyAux
    property = porosity
    variable = porosity
  [../]
  [./permeability]
    type = PorousFlowPropertyAux
    property = permeability
    column = 0
    row = 0
    variable = permeability
  [../]
[]

[Kernels]
  [./precipitation_dissolution]
    type = PorousFlowPreDis
    mineral_density = 1000.0
    stoichiometry = 1
    variable = tracer_concentration
  [../]
[]

[BCs]
  [./constant_injection_of_tracer]
    type = PorousFlowSink
    variable = tracer_concentration
    flux_function = -5E-3
    boundary = injection_area
  [../]
  [./constant_outer_porepressure]
    type = PresetBC
    variable = porepressure
    value = 0
    boundary = rmax
  [../]
[]

[Modules]
  [./FluidProperties]
    [./the_simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 2E9
      viscosity = 1.0E-3
      density0 = 1000.0
    [../]
  [../]
[]

[Materials]
  [./porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
    chemical = true
    initial_mineral_concentrations = initial_and_reference_conc
    reference_chemistry = initial_and_reference_conc
  [../]
  [./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
  [../]
  [./precipitation_dissolution]
    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
  [../]
  [./mineral_concentration]
    type = PorousFlowAqueousPreDisMineral
  [../]
[]

[Preconditioning]
  active = basic
  [./basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  [../]
  [./preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]

[Outputs]
  exodus = true
[]
(modules/porous_flow/examples/tutorial/07.i)

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): (2) 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
  [../]
[]

[Preconditioning]
  active = basic
  [./basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  [../]
  [./preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]

[Outputs]
  exodus = true
[]
(modules/porous_flow/examples/tutorial/07.i)

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

# Darcy flow with a tracer that precipitates causing mineralisation and porosity changes and permeability changes
[Mesh]
  type = AnnularMesh
  dim = 2
  nr = 10
  rmin = 1.0
  rmax = 10
  growth_r = 1.4
  nt = 4
  tmin = 0
  tmax = 1.57079632679
[]

[MeshModifiers]
  [./make3D]
    type = MeshExtruder
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
  [../]
  [./shift_down]
    type = Transform
    transform = TRANSLATE
    vector_value = '0 0 -6'
    depends_on = make3D
  [../]
  [./aquifer]
    type = SubdomainBoundingBox
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    depends_on = shift_down
  [../]
  [./injection_area]
    type = ParsedAddSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    depends_on = 'aquifer'
  [../]
  [./rename]
    type = RenameBlock
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    depends_on = 'injection_area'
  [../]
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [./porepressure]
  [../]
  [./tracer_concentration]
  [../]
[]

[PorousFlowFullySaturated]
  porepressure = porepressure
  coupling_type = Hydro
  gravity = '0 0 0'
  fp = the_simple_fluid
  mass_fraction_vars = tracer_concentration
  number_aqueous_kinetic = 1
  temperature = 283.0
[]

[AuxVariables]
  [./eqm_k]
    initial_condition = 0.1
  [../]
  [./mineral_conc]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./initial_and_reference_conc]
    initial_condition = 0
  [../]
  [./porosity]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./permeability]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[AuxKernels]
  [./mineral_conc]
    type = PorousFlowPropertyAux
    property = mineral_concentration
    mineral_species = 0
    variable = mineral_conc
  [../]
  [./porosity]
    type = PorousFlowPropertyAux
    property = porosity
    variable = porosity
  [../]
  [./permeability]
    type = PorousFlowPropertyAux
    property = permeability
    column = 0
    row = 0
    variable = permeability
  [../]
[]

[Kernels]
  [./precipitation_dissolution]
    type = PorousFlowPreDis
    mineral_density = 1000.0
    stoichiometry = 1
    variable = tracer_concentration
  [../]
[]

[BCs]
  [./constant_injection_of_tracer]
    type = PorousFlowSink
    variable = tracer_concentration
    flux_function = -5E-3
    boundary = injection_area
  [../]
  [./constant_outer_porepressure]
    type = PresetBC
    variable = porepressure
    value = 0
    boundary = rmax
  [../]
[]

[Modules]
  [./FluidProperties]
    [./the_simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 2E9
      viscosity = 1.0E-3
      density0 = 1000.0
    [../]
  [../]
[]

[Materials]
  [./porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
    chemical = true
    initial_mineral_concentrations = initial_and_reference_conc
    reference_chemistry = initial_and_reference_conc
  [../]
  [./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
  [../]
  [./precipitation_dissolution]
    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
  [../]
  [./mineral_concentration]
    type = PorousFlowAqueousPreDisMineral
  [../]
[]

[Preconditioning]
  active = basic
  [./basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  [../]
  [./preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]

[Outputs]
  exodus = true
[]
(modules/porous_flow/examples/tutorial/07.i)

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
  [../]
  [./precipitation_dissolution]
    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
  [../]
  [./mineral_concentration]
    type = PorousFlowAqueousPreDisMineral
  [../]
[]

[Preconditioning]
  active = basic
  [./basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  [../]
  [./preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]

[Outputs]
  exodus = true
[]
(modules/porous_flow/examples/tutorial/07.i)

Instead of using a set of PresetBC 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]
  [./constant_injection_of_tracer]
    type = PorousFlowSink
    variable = tracer_concentration
    flux_function = -5E-3
    boundary = injection_area
  [../]
  [./constant_outer_porepressure]
    type = PresetBC
    variable = porepressure
    value = 0
    boundary = rmax
  [../]
[]
(modules/porous_flow/examples/tutorial/07.i)

Finally, a set of AuxVariables (some employing the handy PorousFlowPropertyAux) is used to record useful information:

[AuxVariables]
  [./eqm_k]
    initial_condition = 0.1
  [../]
  [./mineral_conc]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./initial_and_reference_conc]
    initial_condition = 0
  [../]
  [./porosity]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./permeability]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[AuxKernels]
  [./mineral_conc]
    type = PorousFlowPropertyAux
    property = mineral_concentration
    mineral_species = 0
    variable = mineral_conc
  [../]
  [./porosity]
    type = PorousFlowPropertyAux
    property = porosity
    variable = porosity
  [../]
  [./permeability]
    type = PorousFlowPropertyAux
    property = permeability
    column = 0
    row = 0
    variable = permeability
  [../]
[]
(modules/porous_flow/examples/tutorial/07.i)

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.

Start | Previous | Next