19-aqueous_geochemistry.Rmd

# Aqueous Geochemistry 

## Contributors

Matthew R. Hipsey and S. Ursula Salmon

## Overview

The module allows the user to choose from a wide array of geochemical 'components' (e.g., anions, cations, metals) to simulate, and then solves for the equilibrium speciation of the solution. The module includes an advanced solver, and provides pH and other solution properties, and optional mineral or gas phases. Relevant kinetic transformations (e.g. redox reactions) between simulated components are also included. Metals may be simulated as part of the model and can also be active within the biological cycles. The module is flexible and simple in its configuration and can be used in the water column and in the sediment if the dynamic sediment diagenesis model (see Section ?) is also being simulated. The module dynamically links with the other biogeochemical processes within AED, such that any biological activity (e.g. algal nutrient uptake and photosynthesis, microbial respiration, and nitrification) will dynamically affect the aqueous speciation and pH. Details of the equilibrium and non-equilibrium processes included are described in the adjacent tab.

<center>
```{r 19-pic1, echo=FALSE, fig.cap="Conceptual diagram highlighting key variables and interactions within the model", out.width = '95%'}
knitr::include_graphics("images/19-aqueous_geochemistry/image1.png")
```
</center>

## Model Description

The aqueous geochemistry sub-model was developed to simulate aqueous speciation and solubility equilibrium control of pure-phase minerals using an approach similar to **PHREEQC** and as reported previously in Salmon et al (in review). Each dissolved component, $X$, and mineral phase, $PP$, is a state variable within the model and is mass-conserved. In addition to external loading from inflows, the model includes a configurable dissolved sediment flux term for the dissolved species and a sedimentation term for the particulates. If we remove the hydrodynamic terms (e.g., inflows/outflows, mixing) for brevity, the differential equation is defined for the dissolved state variables as:


###	Process Descriptions


#### Aqueous Speciation & Solubility Equilibrium Control {-}

The aqueous geochemistry model is conceptually similar to other equilbrium codes. The model defines geochemical reactions in terms of components. An aqueous species $C$ is created as a product of reaction between components $A$ and $B$, as shown in the following example:

<center>
\begin{equation}
aA + bB \Longleftrightarrow cC
(\#eq:geochem1)
\end{equation}
</center>

where $a$, $b$ and $c$ are the stoichiometric coefficients. The reaction shown above is able to proceed back and forward depending on the prevailing conditions of the aqueous solution. The equilibrium of the reaction is a function of the standard free energy, and from this the familiar mass-action expression is defined:

<center>
\begin{equation}
K_C=\frac{[C]^c}{[A]^a[B]^b}
(\#eq:geochem2)
\end{equation}
</center>

where $K_{C}$ is denoted the equilibrium constant. In any natural aquatic system, numerous components are present and there are therefore hundreds of reactions similar to Eq. above that will be close to equilibrium. The many reactions form a complex and interdependent set of simultaneous expressions which can be solved numerically, discussed next.

For each component, $X_{x}$, (generally components are equivalent to elements, except for the case of elements with multiple redox states, in which case each redox state is assigned as a unique component), a set of $M$ components is defined, $X \in \{X_1, X_2, ..., X_M\}$. Combination of all the components into an aqueous solution results in numerous reactions taking place of the form of Eq. .., and we define the resultant set of aqueous species that result as: $C \in \{C_1, C_2 ..., C_N\}$, where $N$ is the total number produced.
 
The `activity' of an aqueous species is not equivalent to its molality due to the nonideality of aqueous solutions; a species activity, \tilde{C} is related to its molality according to the activity coefficient, \gamma:

<center>
\begin{equation}
  [C_i] = \tilde{C_i} = \gamma_i C_i
(\#eq:geochem3)
\end{equation}
</center>

Numerous methods exist for estimating a species activity coefficient. The first is known as the *Davies* equation:

<center>
\begin{equation}
  \log \gamma_i = -c_1z_i^2\left(\frac{\sqrt{\mu}}{1+\sqrt{\mu}}-0.3\mu \right)
(\#eq:geochem4)
\end{equation}
</center>

and the other one used in the AED2 solution is the Debye-Huckel equation:

<center>
\begin{equation}
  \log \gamma_i = -\frac{c_1z_i^2\sqrt{\mu}}{1+c_2c_3\sqrt{\mu}} +c_4\mu
(\#eq:geochem5)
\end{equation}
</center>

where $z_{i}$ is the ionic charge of the $i^{\text{th}}$ species $c_{1-4}$ are constants that depend upon properties such as the temperature of the solution. The activity of a species is dependent on the ionic strength, $\mu$, a solution property defined as:

<center>
\begin{equation}
  \mu = \frac{1}{2}\sum_i^N{z_i^2\frac{C_i}{W_{aq}}}
(\#eq:geochem6)
\end{equation}
</center>

where $W_{aq}$ is mass of water in the aqueous phase.

To numerically solve the system of equations, the total number of moles of any component is estimated by adding up the amount contained in each species:

<center>
\begin{equation}
  T_j = \sum_{i=1}^N{a_{ij}C_i} \:\:\:\:\:\:\:\:\:\:\: j=1..M
(\#eq:geochem7)
\end{equation}
</center>

where $))aix$ is the stoichiometric coefficient of component $x$ in species $i$. The activity of each species is calculated using the mass-action equations, which are defined more generically as:

<center>
\begin{equation}
  \tilde{C}_i = K_i \prod_{j=1}^M{\tilde{X}_j^{a_{ij}}} \:\:\:\:\:\:\: j=1..N
(\#eq:geochem8)
\end{equation}
</center>

Taking the logarithm of both sides of Equation above leaves:

<center>
\begin{equation}
  log \tilde{C}_i = \log K_i + \sum_{j=1}^M{a_{ij}\log \tilde{X}_j} \:\:\:\:\: j=1..N
(\#eq:geochem9)
\end{equation}
</center>

which can be defined in matrix notation as:

<center>
\begin{equation}
  \tilde{C}^* = K^*+A\tilde{X}^*
(\#eq:geochem10)
\end{equation}
</center>

where $*$ denotes a vector quantity. The solution is achieved by iterating using a Newton-Rhapson scheme, which aims to minimize the residual, $R$, between the estimated component molalities ($T$) and the known values:

<center>
\begin{equation}
  R_x = \sum_{i=1}^{N}{a_{ix}C_i-T_x}\:\:\:\:\: i=1..N
(\#eq:geochem11)
\end{equation}
</center>

In matrix notation:

<center>
\begin{equation}
  R=A^{-1}C-T
(\#eq:geochem12)
\end{equation}
</center>

The Newton-Rhapson technique iterates towards a solution by employing a derivative (termed Jacobian) matrix, $J$:

\begin{equation}
  R=J \Delta X
(\#eq:geochem13)
\end{equation}

which is solved for $\Delta X$ at each iteration. The new estimates for $X$ are then updated according to:

<center>
\begin{equation}
  X^{n+1} = X^n + \Delta X
(\#eq:geochem14)
\end{equation}
</center>

and the whole scheme proceeds until the solution has sufficiently converged, $R<\epsilon$, where $\epsilon$ is a pre-defined convergence criterion ($10^{-10}$ in this code). The Jacobian matrix is defined as:

<center>
\begin{equation}
  J = \frac{dR_x}{dX_k}
(\#eq:geochem15)
\end{equation}
</center>

which can be expanded to:

<center>
\begin{equation}
  \frac{dR_x}{dX_k} = \sum_{i=1}^{N}{a_{ik}}\frac{dC_i}{dX_k}-\frac{dT_x}{dX_k}
(\#eq:geochem16)
\end{equation}
</center>

using Eq. .... The derivatives are then defined according to:

<center>
\begin{equation}
  \frac{dC_i}{dX_k} = a_{ik}\frac{C_i}{X_k}
(\#eq:geochem17)
\end{equation}
</center>

and

<center>
\begin{equation}
  \frac{dT_x}{dX_k} = 0
(\#eq:geochem18)
\end{equation}
</center>

which are the substituted into Eq. .. to leave

<center>
\begin{equation}
  \frac{dR_x}{dX_k} = \sum_{i=1}^{N}{a_{ik}\:a_{ik}}\frac{C_i}{X_k}.
(\#eq:geochem19)
\end{equation}
</center>

Therefore for a given total number of moles of each element in each computational cell, using this scheme the model solves for the activity of each aqueous species, ionic strength and $pH$. The charge balance variable (denoted $ubalchg$) is also subject to advection and mixing as all other state variables, and importantly, an estimate of $ubalchg$ must be provided at any inflow boundaries. If this is not done, then the charge balance will be compromised, and will manifest in a poor $pH$ prediction.

The Newton-Rhapson scheme described above is implemented by making use of the Simplex numerical solver (Barrodale and Roberts, 1980; Parkhurst and Appelo, 1999). This is because AED2 also allows for simulation of pure phase (i.e. non-aqueous phase) minerals. When minerals can precipitate and dissolve, properties such as the ionic strength ($MU$), activity of water ($XH_{2}O$) and charge balance ($CV$) may be non-conservative in the aqueous solution. These are therefore included as unknowns in the simplex solver, which uses the Newton-Rhapson technique to solve the following matrix:

<center>
\begin{equation}
  \begin{array}{lll}
    \begin{pmatrix}
      R_{PP}\\
      R_{X_{1}}\\
      \vdots\\
      R_{X_{N}}\\
      R_{MU}\\
      R_{H_2O}\\
      R_{CB}\\
      R_{W_{aq}}\\
      X_{PP}\\
      \end{pmatrix}
      =&
      \begin{pmatrix}
      \pd{R_{PP}}{\text{ln}(\tilde{X}_1)} & \dots & \pd{R_{PP}}{\text{ln}(\tilde{X}_N)} & \pd{R_{PP}}{d\mu} & & \\
      \pd{R_{X_1}}{\text{ln}(\tilde{X}_1)} &  &  &  &  & \\
       &  &  &  &  & \\
      \vdots &  &  &  &  & \\
       &  &  &  &  & \\
       &  &  &  &  & \pd{R_{Waq}}{X_{PP}} \\
       &  &  &  &  & 1 \\
       \end{pmatrix}
       &
       \begin{pmatrix}
       d\text{ln}(\tilde{X}_1)\\
       R_{X_{1}}\\
       \vdots\\
       d\text{ln}(\tilde{X}_N)\\
       d\mu \\
       d\text{ln}(\tilde{X}_{H_2O})\\
       d\text{ln}(\tilde{X}_{H^+})\\
       d\text{ln}(W_{aq})\\
       dX_{PP}\\
    \end{pmatrix}
  \end{array}
(\#eq:geochem20)
\end{equation}
</center>

The operation of the simplex solver is complicated and the reader is referred to Barrodale and Roberts (1980) for a detailed account. The aqueous species used in the simulations is based on those from Nordstrom et al. (1990), termed the **WATEQ4F** database.


<!-- #### Redox Transformations {-} -->

<!-- #### Other Geochemistry Sources/sinks (optional) {-} -->




###	Variable Summary

#### State Variables {-}

```{r 19-statevariables, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
state_variables <- read.csv("tables/19-aqueous_geochemistry/state_variables.csv", check.names=FALSE)
kable(state_variables,"html", escape = F, align = "c", caption = "State variables") %>%
  kable_styling(full_width = F, position = "center", font_size =12) %>%
  column_spec(1, width_min = "13em") %>%
  column_spec(2, width_min = "12em") %>%
  column_spec(3, width_min = "12em") %>%
  column_spec(4, width_min = "12em") %>%
  column_spec(5, width_min = "13em") %>%
  row_spec(1:1, background = 'white') %>%
  scroll_box(width = "770px", 
             fixed_thead = FALSE) 
```

###	Parameter Summary
```{r 19-parameterconfig, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
parameter_config <- read.csv("tables/19-aqueous_geochemistry/parameter_config.csv", check.names=FALSE)
kable(parameter_config,"html", escape = F, align = "c", caption = "Parameters and configuration",
      bootstrap_options = "hover") %>%
  kable_styling(parameter_config, bootstrap_options = "hover", 
                full_width = F, position = "center",
                font_size = 12) %>%
  column_spec(1, width_min = "6em") %>%
  column_spec(2, width_min = "6em") %>%
  column_spec(3, width_min = "6em") %>%
  column_spec(4, width_min = "6em") %>%
  column_spec(5, width_min = "6em") %>%
  column_spec(6, width_min = "6em") %>%
  column_spec(7, width_min = "6em") %>%
  row_spec(1:19, background = 'white') %>%
  scroll_box(width = "770px", height = "550px",
             fixed_thead = FALSE) 
```
<br>

<!-- ###	Optional Module Links -->
<!-- ###	Feedbacks to the Host Model -->


## Setup & Configuration

### Setup Example

An example `nml` block for the geochemistry module is shown below:

```{fortran, eval = FALSE}
&aed_geochemistry
   geochem_file   = './AED2/aed_geochem_pars.dat'
   !-- dissolved components
     num_components = 10
     dis_components = 'DIC', 'FeII', 'FeIII' ,'Al', 'Cl', 'SO4', 'Na', 'K', 'Mg', 'Ca','H2S'
     component_link = 'CAR_dic','','','','','','','','','','','',''
     Fsed_gch       = 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.
   !-- mineral components
     num_minerals   = 3
     the_minerals    = 'Gibbsite','FeOH3A','Calcite'
     mineral_link    =  '',   '',  ''
     w_gch           = 0.0, -1.0, 0.0
   !-- iron redox
     Riron_red = 0.0
     Kiron_red = 100.
     theta_iron_red = 1.07
     Riron_aox = 0.0  !1 = 100% per day decay
     Riron_box = 0.5
     theta_iron_ox = 1.06
   !-- advanced
     simEq = .true.
   ! dis_initial    = 0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1
   ! speciesOutput  = 'CO2'
   ! min_initial    = 0.01,0.01,0.01,0.01
  /
```



<br>


```{block2, gch-text, type='rmdnote'}
In addition to adding the above code block to `aed.nml`, users must also supply a valid AED geochemistry reaction database file (`aed_geochem_pars`). This file must be supplied in the specified `DAT` format. 
```

<br>


<!-- ## Case Studies & Examples -->
<!-- ###	Case Study -->
<!-- ###	Publications -->