Geochemical Transport#

Geochemical transport capabilities are implemented in a few different PKs.

Transport#

src/physics/ats/src/pks/transport/transport_ats.hh

This PK solves for transport of chemical species in water. It may optionally be paired with a PK for reactions, typically in an operator-split form, to solve reactive transport.

The advection-dispersion equation for component i may be written as:

\[\frac{\partial (\Theta \xi_i)}{\partial t} = - \boldsymbol{\nabla} \cdot (\boldsymbol{q} \xi_i) + \boldsymbol{\nabla} \cdot (\Theta \, (\boldsymbol{D^*}_l + \tau \boldsymbol{D}_i) \boldsymbol{\nabla} \xi_i) + Q_s,\]

The primary variable for this PK is \(\xi_i\), the mole ratio of a specie \(i\) in the aqueous phase, with units of [mol i / mol H2O].

This seems a bit odd to most people who are used to reactive transport codes where the primary variable is total component concentration \(C\), in units of [mol i / L]. The reason for this differences is to better enable variable-density problems. Note that the two are related:

\[C_i = n_l \xi_i\]

for molar density of liquid water \(n_l\).

For reactive transport problems, both concentration and mole ratio are output in the visualization file. For transport-only problems, concentration may be output by adding a total component concentration evaluator that multiplies the two quantities using an “evaluator type” = “multiplicative evaluator”.

“PK type” = “transport ATS”

pk-transport-ats-spec

  • “domain name” [string] “domain” specifies mesh name that defines the domain of this PK.

  • “component names” [Array(string)] optional No default. Provides the names of the components that will be transported. Note that these must be provided by the user if transport is used alone, or are provided by the geochemical engine if reactions are to be used.

  • “number of aqueous components” [int] -1 The total number of aqueous components. Default value is the length of “component names”

  • “boundary conditions” [transport-bc-spec] Boundary conditions for transport are dependent on the boundary conditions for flow. See Transport-specific Boundary Conditions

  • “molecular diffusivity [m^2 s^-1]” [molecular-diffusivity-spec] See below.

  • “tortuosity [-]” [tortuosity-spec] See below.

  • “source terms” [transport-source-spec-list] Provides solute source.

Math and solver algorithm options:

  • “diffusion” [pde-diffusion-typedinline-spec] Diffusion drives the distribution. Typically we use finite volume here, but mimetic schemes may be used. See Diffusion

  • “diffusion preconditioner” [pde-diffusion-typedinline-spec] Inverse of the above. Likely only Jacobian term options are needed here, as the others default to the same as the “diffusion” list. See Diffusion.

  • “inverse” [inverse-typed-spec] Inverse method for the diffusion-dispersion solve. See Inverse.

  • “cfl” [double] 1.0 Courant-Friedrichs-Lewy condition, a limiter on the timestep size relative to spatial size. Must be <= 1.

  • “advection spatial discretization order” [int] 1 Defines accuracy of the spatial discretization in the advection term. It permits values 1 or 2. Default value is 1 (donor upwind) but 2 (a limiter scheme) is much less numerically diffusive, and recommended for most cases.

  • “temporal discretization order” [int] 1 Defines accuracy of temporal discretization. It permits values 1 (forward Euler) and 2 (a Runga-Kutta scheme).

  • “reconstruction” [reconstruction-spec] collects reconstruction parameters for use in reconstructing the velocity field for 2nd order advection schemes. See Reconstructions.

KEYS

  • “primary variable” “molar_ratio” [mol C mol H2O^-1]

  • “liquid water content” “water_content” This variable is a multiplier in in the accumulation term. This is often just “water_content”, but not in frozen cases, where it must be defined by the user (typically as \(\phi s n_l |V|\) in the subsurface, or \(h n_l |V|\) on the surface.

  • “molar density liquid” [mol H2O m^-3]

  • “water flux” The face-based water flux in [mol H2O s^-1].

  • “water source” Defines the water injection rate [mol H2O m^-2 s^-1] in surface and [mol H2O m^-3 s^-1] in subsurface) which applies to concentrations specified by the “geochemical conditions”. Note that if this PK is coupled to a surface flow PK, the unit of the water source there must be in [mol m^-2 s^-1], not in [m s^-1] as is an option for that PK (e.g. “water source in meters” must be set to “false” in the overland flow PK).

    The injection rate of a solute [molC s^-1], when given as the product of a concentration and a water source, is evaluated as:

    Concentration [mol C L^-1] \(\times\) 1000 [L m^-3] of water \(\times\) water source [mol H2O m^-3 s^-1] \(\times\) volume of injection domain [m^3] \(/\) molar density of water [mol H2O m^-3]

Note, this is not dispersion, but strictly isotropic diffusion.

molecular-diffusivity-spec

For each aqueous component, a single value is provided for molecular diffusivity, e.g.

“COMPONENT_NAME” [double] value [m^2 s^-1]

tortuosity-spec

  • “aqueous” [double] Defines tortuosity for calculating diffusivity of liquid solutes, [-].

  • “gaseous” [double] Defines tortuosity for calculating diffusivity of gas solutes, [-]. Not currently implemented!

transport-source-spec

  • “component mass source” [list] Defines solute source injection rate.

    • “spatial distribution method” [string] One of:

      • “volume”, source is considered as extensive quantity [molC s^-1] and is evenly distributed across the region.

      • “none”, source is considered as intensive quantity. [molC m^-2 s^-1] in surface and [molC m^-3 s^-1] in subsurface

    • “geochemical” [list] Defines a source by setting solute concentration for all components (in moles/L) and an injection rate given by the water source. Currently, this option is only available for Alquimia provided geochemical conditions.

      • “geochemical conditions” [Array(string)] List of geochemical constraints providing concentration for solute injection.

Chemistry#

Reactive transport is done by coupling a Transport PK with chemistry using some variation on a Reactive Transport MPC. Currently, ATS only implements chemistry through Alquimia, which calls out to either PFloTran or Crunch for reaction support.

src/pks/chemistry/Alquimia_PK.hh

The Alquimia PK leverages the Alquimia interface to implement chemical reactions by calling out to a chemical engine.

This PK, when not coupled, solves the following set of ODEs:

\[\frac{\partial C_ij}{\partial t} = \sum_{k=1..nreactions} R_ik(C_{i=1..ncomponents,j})\]

for \(k\) indexing reactions, \(i\) indexing components, and \(j\) indexing grid cells. The total component concentrations \(C\) are in [kg L^-1].

The details of how to parameterize and specify a reaction network are specific to each geochemical engine, of which currently only PFloTran and Crunch are supported. Corresponding input files for these engines must be provided. Examples are available in Amanzi’s user guide tests and ATS’s regression tests.

“PK type” = “chemistry alquimia”

pk-chemistry-alquimia-spec

  • “timestep controller” [timestep-controller-typed-spec] Timestep control parameters.

  • “engine” [string] One of:

    • “CrunchFlow”

    • “PfloTran”

  • “engine input file” [string] Path to the engine’s input file.

KEYS:

  • “total component concentration”

  • “mass density liquid”

  • “porosity”

  • “saturation liquid”

  • “temperature”

  • “mineral volume fractions”

  • “mineral specific surface area”

  • “mineral rate constant”

  • “surface site density”

  • “total_sorbed”

  • “isotherm kd”

  • “isotherm freundlich n”

  • “isotherm langmuir b”

  • “first order decay ratae constant”

  • “cation exchange capacity”

  • “aux data”

  • “pH”

See the Alquimia documentation for the needed input spec for a given engine.

<ParameterList>  <!-- parent list -->
<ParameterList name="alquimia PK">
  <Parameter name="engine" type="string" value="PFloTran"/>
  <Parameter name="engine input file" type="string" value="tritium.in"/>
  <ParameterList name="timestep controller">
  <Parameter name="timestep control type" type="string" value="fixed"/>
  <ParameterList name="fixed parameters">
    <Parameter name="timestep [s]" type="double" value="1.0e-02"/>
  </ParameterList>
</ParameterList>

Salinity Transport#

Salinity transport is not itself a PK, but may be implemented as a variation on Transport by supplying equations of state that are a function of the concentration of salt. Examples of this are provided in the regression tests and demos.

Sediment Transport#

src/physics/ats/src/pks/transport/sediment_transport_pk.hh

Document me @dasvyat