23-sediment_biogeochemistry.Rmd

# (PART) AED Benthic Modules {.unnumbered}

# Sediment Biogeochemistry

## Contributors

Daniel Paraska, Matthew R. Hipsey

## Overview

This module is a sediment reactive transport model, based on a 1D approximation of the sediment solid and pore-water profiles. Each active sediment zone (or column) is discretized into a user defineable number of layers that start at thicknesses of around 1 mm at the sediment-water interface and increase exponentially down to a pre-defined sediment depth. The model resolves in each layer both physical processes (e.g. pore-water diffusion or bioturbation) and chemical processes (e.g. redox transformations).

<center>

```{r SedOverview-1, echo=FALSE, fig.cap="Schematic of the main physical and chemical processes that cause chemical concentration and flux change in the sediment and across the sediment-water interface, and therefore are included in most sediment diagenesis models.", out.width = '85%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Boxes&Arrows/Overview-05.png")
```

</center>

Under some conditions, the sediment stores can release nutrients to the water column, while under other conditions, the sediment can remove nutrients, through burial over the long-term, and through processes such as denitrification in the surface layers. The fine balance that controls the conditions under which the sediment will store, release or remove nutrients is largely governed by the aerobic state of the sediment pore water, and the amount of reactive organic matter fuelling the reactions. The depth-resolved sediment model accounts for mixing from the hydrodynamic model into the upper sediment layers, and then calculates whether organic matter is consumed aerobically, through denitrification or, deeper down, through sulfate reduction or even methanogenesis.

Sediment early-diagenesis models are complex environmental reactive transport simulation tools. The meta-analysis by Paraska et al. (2014) discussed the history of their evolution to these complex configurations, in which the original models of Boudreau (1996), Van Cappellen and Wang (1996) and Soetaert et al. (1996) were taken up and applied in many contexts by new modellers, who added new features and extended their capabilities, or discarded old features as required. The meta-analysis also identified the major challenges associated with developing new sediment diagenesis models. Here, the AED modelling package for sediment biogeochemistry is presented, CANDI-AED, which is an extension of the Approach 1 models, but reengineered and augmented with new model approaches and capabilities as a way to address some of these challenges.

Paraska et al. (2015) outlined the significance and uncertainty associated with different parameterisation approaches of organic matter dynamics. In these cases, simulations were run to test the significance of different theoretical approaches and model structural assumptions, using an idealised model setup with only primary oxidation reactions and no physical processes or spatial resolution. The true impact of these different model approaches within a spatially-resolved model, accounting for all of the advection, diffusion and secondary reaction processes, however, is yet to be determined and it is unclear whether some formulations may suit some application contexts better than others. Therefore there is a need for a fully flexible model structure that can include these different organic matter breakdown parameterisations and allow users to assess critically the alternative approaches. In addition, other aspects related to secondary redox reactions, mineral reactions, precipitation and adsorption should similarly be subject to comparative assessments.

The model included in AED aimed to address challenges of building a generic and full-featured, open-source model code with the flexibility to do the following:

-   set different kinetic rate equation approaches
-   set different organic matter pools and breakdown processes
-   use standard inhibition or thermodynamic limits on primary oxidation
-   optionally use manganese, iron and iron sulphide reactions
-   simulate adsorbed metals and nutrients
-   simulate calcium, iron and manganese carbonates
-   connect the boundary to either another model, a programmed file or a fixed concentration

Therefore the numerical model presented in this module has many optional features and alternative parameterisations for key processes, without mandating their inclusion in the calculations or enforcing a fixed model structure.

The sediment model CANDI-AED presented here is implemented as an optionally configurable module in the AED model library. Through the model coupling approach it may be applied with any of the hydrodynamic models linked to AED (e.g. GLM, Tuflow), or alternatively, options to run in isolation are also possible. This chapter provides a scientific description of the model and describes attributes of the model associated with its practical implementation and operation. A case-study of the model framework is also demonstrated.

The major output of CANDI-AED includes concentration depth profiles and fluxes across the sediment-water interface (figure \@ref(fig:SedOverview-2)).

<center>

```{r SedOverview-2, echo=FALSE, fig.cap="Examples of the final outputs after running CANDI-AED: concentration-depth profiles and flux-time plots.", out.width = '100%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Outputs/Overview-17.png")
```

</center>

## Model Description

### General description

The heart of this model is the reaction, diffusion, advection model of Berner (1980), which was implemented as the Carbon and Nutrient Diagenesis model of Boudreau (1996), later further developed as the C.CANDI code (Luff et al. 2000). The CANDI-AED implementation, however, has evolved from the original code, and includes extensions related to the treatment of organic matter, the simulation of the geochemical conditions known to influence the diagenetic equations, extensions for nutrients and trace metals, and dynamics at the sediment-water interface. However, the core organic matter breakdown equations (and their numerical solution) remains similar as the original descriptions presented in Boudreau (1996), and to other similar sediment models. An overview of the model is shown in Figure \@ref(fig:23-pic2).

<center>

```{r 23-pic2, echo=FALSE, fig.cap="Overview of the chemical processes in CANDI-AED: organic matter transformation and oxidation, and reduction/oxidation, crystallisation, adsorption and precipitation reactions of inorganic by-products. Most of the processes are triggered by the input of POM at the sediment-water interface.", out.width = '85%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Boxes&Arrows/FullCandi-02.png")
```

</center>

<br>

The model is based on the advection-dispersion reaction equation for the concentration of dissolved and particulate substances. For dissolved substances $C_d$, the balance equation is defined as:

<center>

```{=tex}
\begin{equation}
\frac{\phi \delta C_d}{\delta t} = \overbrace{ D_{B}\frac{\delta^{2} C_d}{\delta x^{2}} + \phi D_{S}\frac{\delta^{2} C_d}{\delta x^{2}} }^\text{biodiffusion and molecular diffusion} - \underbrace{u \frac{\delta C_d}{\delta x}}_\text{advection (flow)} + \overbrace{\alpha(C_{d_0}-C_d)}^\text{irrigation}  + \underbrace{\phi \sum R_d}_\text{reactions} + \color{brown}{S}
(\#eq:sdg-1)
\end{equation}
```
</center>

where the left hand side (LHS) is the "unsteady" term ($C_d$ change in time), the first term on the right hand side (RHS) is the dispersion/mixing term, the second term on the RHS is the advection/movement term, and $R_d$ denotes a generic reaction term. An optional $S$ term is included to represent sources from other modules (e.g., seagrass root injection of $O_2$).

For particulate (solid) substances, $C_s$:

<center>

```{=tex}
\begin{equation}
\frac{(1-\phi)\rho \delta C_s}{\delta t} = \overbrace{ D_{B}\frac{\delta^{2}[(1-\phi)\rho C_s] }{\delta x^{2}} }^\text{biodiffusion} - \underbrace{\omega \frac{\delta [(1-\phi)\rho C_s]}{\delta x}}_\text{advection (sedimentation)}  + \overbrace{(1-\phi) \sum R_s}^\text{reactions} + \color{brown}{S}
(\#eq:sdg-2)
\end{equation}
```
</center>

where $(1-\phi)$ denotes the solid fraction of the sediment, and $R_s$ is a generic reaction term.

The above equations are solved numerically for the simulated set of constituents. The user can define the variables that are included in the $\mathbf{SDG}$ module, such that $C_d$ and $C_s$, are the set of dissolved and particulate variables selected for simulation, respectively. Eight of the variables (compulsory variables) must be requested by the user in the variable setup file, or the model will not run (see the [variables summary](#VariableSummary) below). Three other compulsory variables are always simulated and the user does not need to request them. The other variables can be requested if the user desires them. A number of options are available for resolving the physical processes, including the rate of diffusion, advection, irrigation and the boundary condition options.

In addition to physical processes, the CANDI-AED model considers two types of chemical reactions - the slow, kinetically controlled reactions, and the fast thermodynamically based equilibrium reactions. The latter are simulated in the sediment through appropriate configuration of the [geochemistry reactions](#EquilibriumGeochemistry); the configuration of the equilibrium model will apply to both the water and each of the sediment layers. The kinetically controlled reactions are mostly microbially-mediated and include the reactions for [organic matter breakdown](#SedPrimaryRedoxReactions) and eventual oxidation, the [re-oxidation of various by-products](#SedSecondaryRedoxReactions) and the dynamics of the metal sulfides. These reactions can be complex and are outlined in further detail in the next sections.

### Process Descriptions

#### Sediment model domain

The sediment model is discretised into a user-definable number of depth layers (`maxnpts`) down to a pre-defined sediment depth(`xl`). The grid of layers is set to have either even spacing (`job` = 1) or to be exponentially increasing (`job` = 0), where the layers sizes have a thickness of a few mm at the sediment-water interface and which increase exponentially down into the sediment (figure \@ref(fig:Init-1)). When the spacing increases exponentially, the first two layers are hard-coded to be a miniumum of 0.25 cm, in order to avoid numerical instability from sharp concentration differences over distances that are too small.

This fixed depth of sediment has a concentration of all variables at all depth layers, and boundary conditions at the upper and lower ends of the domain.

```{r Init-1, echo=FALSE, fig.cap="Initialisation of the depth layers. The number of layers is set by `maxnpts` and the depth of the simulation by `xl`. The setup can have even spacing (left) using parameter `job` = 1 or increasing spacing (right) `job` = 2.",fig.show='hold',fig.align='center' ,out.width = '60%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Initialisation/job1job2-01.png")
```

#### Physical Transport

##### *Solids* {.unnumbered}

CANDI-AED adopts the approach of Boudreau (1996) to advection and dispersion, which is similar to most other diagenesis models. Advection of the solid sediment matrix is occurs at the rate of sediment deposition ($\omega_{00}$ in cm y^-1^). The sediment does not move, however, since the height of the modelling domain is fixed, as more sediment accumulates at the surface, a sediment particle moves downwards. An alternative way to view this is that the modelling domain moves upwards.

```{r Phys1, echo = FALSE, out.width='60%', class = "text-image",fig.show='hold',fig.align='center', fig.cap = "Upon sedimentation, the frame of reference shifts upwards. The bottom area is effectively lost from the modelling domain over time. Particles and porewater are advected downward relative to the modelling domain."}

knitr::include_graphics("images/23-sediment_biogeochemistry/AdvectionDispersion/AdvectionDispersion-05.png")

```

##### *Porewater* {.unnumbered}

For the porewater components, diffusion coefficients are used that are based on free-solution molecular diffusion constants corrected for sediment tortuosity, $θ$. Porewater moves downwards at the same rate as the solids ($\omega_{00}$). A further porewater advection term (`poreflux`) is available, which could represent, for example, pressure on the sediment from a groundwater source. If `poreflux` is positive, then porewater moves upwards, relative to particles. Conversely, if `poreflux` is negative, then the porewater moves downwards relative to the particles. In most simulations, `poreflux` is zero and advection is the transport process. A dynamic `poreflux` can be assigned using the column headings `w00h00` in [swibc.dat](#DynamicOptionsSwibcFile) or `pf` in [deepbc.dat](#DynamicOptionsDeepbcFile) (see setup section below).

```{r Phys2, echo = FALSE, out.width='100%', class = "text-image",fig.show='hold',fig.align='center', fig.cap = "Porewater is advected at the same rate as burial. If `poreflux` is non-zero, an additional rate of advection is applied to porewater. The direction is upwards when `poreflux` is positive, downwards when *poreflux* is negative."}

knitr::include_graphics("images/23-sediment_biogeochemistry/AdvectionDispersion/Poreflux-04.png")

```

Diffusion of solutes is calculated using a diffusivity constant. One default value, which is temperature-dependent, is used for most dissolved variables (table \@ref(tab:Diff-table)). This value is based on the diffusivity of chloride. Twelve other variables have their own values, which are also temperature-dependent. Some also use the function *a*, which can also include salinity and pressure (table \@ref(tab:Diff-table)).

<center>

```{r Diff-table, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
Diff_tab <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="Diffcoef2")
  kable(Diff_tab[,1:2],"html", escape = F, align = "c"
        , caption = "Temperature-dependent calculations for diffusivity constants. All state variables in CANDI-AED use the default value, unless specified in this table. The function *a* is a further correction for salinity and pressure.",bootstrap_options = "hover")%>%
   kable_styling(Diff_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>% 
    
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "20em" ,color="black",bold = F) %>%
  column_spec(2, width_min = "20em" ,color="black") %>%
  # column_spec(3, width_min = "15em" ,color="black") %>%
  # row_spec(1:2, background = 'white')  %>%
 scroll_box(width = "40em", height = "20em",fixed_thead = FALSE)
```

</center>

##### *Porosity* {.unnumbered}

Porosity ($\phi$) is defined as the amount of water per total volume of space, as a real number between 0 and 1. A value for porosity is set at each depth layer during initialisation (see [initialisation section](#SedInitialisation) below) and remains constant throughout the simulation (Figure \@ref(fig:Phys-6)). The porosity array is used in the model to calculate other variables, such as

-   advection
-   poreflux
-   porewater volume
-   solid volume
-   diffusive velocity
-   solid fraction

Porosity is also used in the model to calculate $ps$, $psp$ and $pps$, which are used to convert between solid and dissolved state variables, using these conversions:

```{=tex}
\begin{eqnarray}
     ps  = 1 - \phi
     \\
     (\#eq:poros-1)
     \\
     psp = \frac {ps} {\phi}
     \\
     (\#eq:poros-2)
     \\
     pps = \frac {\phi} {ps}
     (\#eq:poros-3)
\end{eqnarray}
```
Porosity is configured by setting three parameters: the porosity at the sediment-water interface (`p0`), the porosity at the depth where compaction causes constant porosity with depth (`p00`), and the attenuation coefficient on the exponential curve (`bp`) (equation \@ref(eq:init-1)). These parameters are set separately for each zone in `aed_candi_params.csv`. Porosity is the amount of water per total volume of space, as a real number between 0 and 1.

```{=tex}
\begin{eqnarray}
    \phi = (\rho _0 - \rho _{00}) \times e^{-bp \times depth}  + \rho _{00}
(\#eq:init-1)
\end{eqnarray}
```
<br>

```{r Phys-6, echo=FALSE, fig.cap="Example depth profile of porosity.",fig.show='hold',fig.align='center' ,out.width = '60%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Initialisation/Porosity-01.png")
```

\

Porosity and bioirrigation profiles are output in the file `Depths.sed`, along with the depths (in cm), and so the user can examine the shapes of the profiles with depth.

#### Primary Redox Reactions {#SedPrimaryRedoxReactions}

The key chemical process that causes ongoing change in the sediment is the breakdown of organic matter. Organic matter supplies fuel, above all, reduced carbon, to microorganisms living in the sediment.

##### *Organic matter pools and reactivity* {.unnumbered}

Organic matter ($OM$) degradation pathways can include labile and refractory components, and the breakdown pathways simulated are summarized conceptually in Figure \@ref(fig:OMModels). Reactions included in the kinetic component include the hydrolysis of the complex (e.g., high molecular weight) $OM$ pools ($POM_{VR}$, $POM_R$, $DOM_R$, $POM_L$) and transformation of low molecular weight (LMW) $DOM_L$ by oxidants ($O_2$, $NO_3^-$, $MnO_2$, $Fe(III)$ and $SO_4^{2-}$ - the so-called ['terminal metabolism'](#SedTheRedoxSequence) pathway), and the release of resulting nutrients ($NH_4^+$, $PO_4^{3-}$) and reduced by-products ($Mn^{2+}$, $Fe^{2+}$, $N_2$, $H_2S$, $CH_4$) and $CO_2$. Oxidants, nutrients and by-products are all capable of interacting, for example, through re-oxidation of reduced species (outlined in the [secondary redox section](#SedSecondaryRedoxReactions) below).

The user can decide how complex or simple the organic matter breakdown pathway should be, with three options of varying complexity for parameterising the pathways included (Figure \@ref(fig:OMModels)). The first option (`OMModel` = 1) is a common multi-G model in which the $POM$ phases are decomposed straight to $CO_2$ and other breakdown products. Here $POM$ is a variable that is not precisely defined, and its components (such as C, N and P) are assigned by parameters based on a user-defined stoichiometry. `OMModel` 1 has a switch `VCW` (logical) that replaces the kinetic rate constants for $POM$ with a depth-dependent breakdown rate, as used by Van Cappellen and Wang (1996).

The second option (`OMModel` = 2) is another '2G' model with both particulate and dissolved organic matter ($POM$ and $DOM$) phases included and parameterisation hydrolysis of $POM$ to $DOM$, and then $DOM$ to $CO_2$ and other breakdown products. The $POM$ and $DOM$ phases consist of three variables each, which trace the reaction and transport of carbon, nitrogen and phosphorus, thereby allowing for variable stoichiometry of organic matter to occur temporally and spatially.

The third option (`OMModel` = 3) has many $POM$ phases, which are all hydrolysed to $DOM$, which then undergoes fermentation and terminal metabolism. This allows the carbon, nitrogen and phosphorus to be calculated precisely before and after a model run, and allows the free energies of the reaction of each phase to be included. This third option is the most detailed and mechanistic, and allows for expansion of more detailed reaction mechanisms to be included, but is recommended only for experienced users.

CANDI-AED also contains options for very unreactive organic matter parameter. The model can decrease the reactivity of the particulate refractory phase towards zero as the concentration approaches a minimum. For example, for $POC_R$, the parameter `pocu` sets the minimum concentration and `KOMPres_C` sets the sensitivity of the decrease, as in equation \@ref(eq:pomr-1) below. The equivalent parameters for N an P are `ponu` and `popu`, and `KOMPres_N` and `KOMPres_P`.

```{=tex}
\begin{eqnarray}

komlim = \frac {PON_R - ponu} {PON_R - ponu + KOMPres_N}
\\
\\
R_{PONR} = (ponr2donr) (komlim)
(\#eq:pomr-1)

\end{eqnarray}
```
(The equation above is a simplified version, where the full equation had checks to prevent the rate going below zero, plus the respiration processes of $MPB$.)

```{r OMModels, echo=FALSE, fig.cap="Three options for different levels of complexity in organic matter breakdown, by setting the `OMModel` switch. Top left – Model in which $POM$ breaks down directly to $CO_2$ and other waste products. Top right - Model with the breakdown rate set by depth. Bottom left – Model in which POM is first hydrolysed to $DOM$ and then oxidised to $CO_2$. Bottom right – Model in which $POM$ is hydrolysed to $DOM$, which can then be fermented and oxidised.", fig.show="hold", out.width="75%" ,fig.align='center'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Boxes&Arrows/OMModel-06.png")
```

<br>

<center>

```{r OMModel-table, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
OMModel_tab <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="OMModel")
# View(sediment_VarsTable)
OMModel_Groups <- unique(OMModel_tab$Group)
  kable(OMModel_tab[,1:3],"html", escape = F, align = "c"
        , caption = "Organic matter breakdown processes. The index *i* in *∑ R~Oxi~* refers to the sequence of *TEA* in table [](@ref(tab:OMChem-table)). For `OMModel` 2, $POM$ and $DOM$ are each three state variables of $POC$, $PON$ and $POP$, and $DOC$, $DON$ and $DOP$.",bootstrap_options = "hover")%>%
  
    
  pack_rows(OMModel_Groups[1],
            min(which(OMModel_tab$Group == OMModel_Groups[1])),
            max(which(OMModel_tab$Group == OMModel_Groups[1])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(OMModel_Groups[2],
            min(which(OMModel_tab$Group == OMModel_Groups[2])),
            max(which(OMModel_tab$Group == OMModel_Groups[2])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(OMModel_Groups[3],
            min(which(OMModel_tab$Group == OMModel_Groups[3])),
            max(which(OMModel_tab$Group == OMModel_Groups[3])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(OMModel_Groups[4],
            min(which(OMModel_tab$Group == OMModel_Groups[4])),
            max(which(OMModel_tab$Group == OMModel_Groups[4])),
            background = '#ebebeb'
            ,color="black") %>%

   kable_styling(OMModel_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 11) %>% 
    
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "18em" ,color="black",bold = F) %>%
  column_spec(2, width_min = "18em" ,color="black") %>%
  column_spec(3, width_min = "18em" ,color="black") %>%
  # column_spec(4, width_min = "15em") %>%
  row_spec(1:4, background = 'white')  %>%
 scroll_box(width = "770px", height = "28em",fixed_thead = FALSE)
```

</center>

##### *C:N:P partitioning* {.unnumbered}

The C:N:P stoichiometry is set differently for each `OMModel`.

Using `OMModel` = 1, the ratios are set separately for the labile and refractory phases, using the *x*, *y* and *z* ratios outlined in table \@ref(tab:OMChem-table). For example, `xlab` is the carbon fraction of $POM_L$, which is usually set to around 106, while `zlab` is the $POM_L$ phosphorus fraction, which is usually set to around 1. These ratios are fixed throughout the simulation. The *x*, *y* and *z* coefficients are used to set the proportions of the other species in the redox equation (such as oxidants and reduction products), the acid-base products (such as $HCO_3^-$ and $CO_2$) and the amounts of free $NH_4^+$ and $PO_4^{3-}$ released from organic matter degradation oxidation.

The *x*, *y* and *z* ratios in table \@ref(tab:OMChem-table) are written in the standard form of the equations in the sediment model literature. The carbon in the organic matter is assumed to be carbon (0) and each carbon atom requires four electrons from the oxidant to form the carbon (IV) in $CO_2$, and therefore *x* carbons require 4*x* electrons. In practice, the model code performs reactions per carbon, and so all of these coefficients are divided by *x*, for all reactants and products. This is done in order to make the reactions `OMModel` 1 equivalent to `OMModel` 2.

Using `OMModel` = 2, the C:N:P ratio is set by the concentrations of separate organic species, such as $POC_L$, $PON_L$ and $POP_L$, as opposed to the $POM$ phases of `OMModel` 1. The ratio is fixed by the boundary and initial conditions, however, if the boundary is dynamic, then the ratio can change through the course of the simulation. Only $DOC_L$ participates in the redox reactions, and so the equivalent of *CH~2~O* in table \@ref(tab:OMChem-table) is $DOC_L$. The release of organic N as free $NH_4^+^$ is calculated $PON_L$ oxidation at `ponl2dic` and the release of organic P as free $PO_4^{3-}$ as $POP_L$ at `popl2dic`.

Using `OMModel` = 3 is similar to `OMModel` = 1. The ratios are applied to $POM_1$, $POM_2$ etc. as `xPOM1`, `yPOM2` etc.

##### *The redox sequence* {#SedTheRedoxSequence .unnumbered}

The terminal redox reaction pathways are the six major pathways that are available in most diagenesis models, and are driven by different organic matter pools, depending on the `OMModel` configuration chosen from the above options. CANDI-AED allows the use of Approach 1 or 2 organic matter oxidation rate equations, as examined in detail in Paraska et al. (2015), and also explained below. The six major pathways are:

-   Aerobic
-   Denitrifying
-   Manganese reduction
-   Iron reduction
-   Sulfate reduction
-   Methanogenesis

The decay of the complex $OM$ types to the LMW $DOM$ required for the heterotrophic bacteria to utilise are all modelled with a simple first-order decay rate. The subsequent reactions for terminal metabolism that describe the breakdown of $OM$ to $CO_2$ and other breakdown products are described in table \@ref(tab:OMChem-table).

CANDI-AED uses a more detailed set of reactions for the denitrifying process, which are described in the [Nitrogen cycling](#SedNitrogenCycling) section below. The chemical equations may be written as in table \@ref(tab:OMChem-table) and the reaction rates for each of these are calculated dynamically based on Monod expressions which mediate the reaction rate according to the concentration of potential oxidants higher in the redox sequence, and the concentration of the available oxidant.\@ref(fig:dev-Grid1)

<br>

<center>

```{r OMChem-table, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
OMChem_tab <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="OMChem3")
# View(sediment_VarsTable)
OMChem_Groups <- unique(OMChem_tab$Group)
  kable(OMChem_tab[,1:3],"html", escape = F, align = "c"
        , caption = "Primary terminal redox reactions. x, y and z are stoichiometric coefficients.",bootstrap_options = "hover")%>%
  
  #   
  # pack_rows(OMChem_Groups[1],
  #           min(which(OMModel_tab$Group == OMModel_Groups[1])),
  #           max(which(OMModel_tab$Group == OMModel_Groups[1])),
  #           background = '#ebebeb'
  #           ,color="black") %>%
  # pack_rows(OMChem_Groups[2],
  #           min(which(OMModel_tab$Group == OMModel_Groups[2])),
  #           max(which(OMModel_tab$Group == OMModel_Groups[2])),
  #           background = '#ebebeb'
  #           ,color="black") %>%
  # pack_rows(OMChem_Groups[3],
  #           min(which(OMModel_tab$Group == OMModel_Groups[3])),
  #           max(which(OMModel_tab$Group == OMModel_Groups[3])),
  #           background = '#ebebeb'
  #           ,color="black") %>%
  # pack_rows(OMChem_Groups[4],
  #           min(which(OMModel_tab$Group == OMModel_Groups[4])),
  #           max(which(OMModel_tab$Group == OMModel_Groups[4])),
  #           background = '#ebebeb'
  #           ,color="black") %>%

   kable_styling(OMChem_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 11) %>% 
    
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "10em" ,color="black",bold = T) %>%
  column_spec(2, width_min = "2em" ,color="black") %>%
  column_spec(3, width_min = "15em" ,color="black") %>%
  # column_spec(4, width_min = "15em") %>%
  row_spec(1:4, background = 'white')  %>%
 scroll_box(width = "770px", height = "40em",fixed_thead = FALSE)
```

</center>

<br>

<br>

<!-- \begin{eqnarray} -->

<!-- \\ -->

<!-- N2O breakdown -->

<!-- \\ -->

<!-- (\#eq:primary-0) -->

<!-- \\ -->

<!-- \overbrace{(CH_2O)_{106}(NH_3)_{16}(H_3PO_4)}^{\textrm{organic matter}} &+& 138O_2 \rightarrow \nonumber -->

<!-- \\ -->

<!-- && 106CO_2 + 16HNO_3 + H_3PO_4 + 122H_2O -->

<!-- \\ -->

<!-- \textrm{Free energy, } \Delta G_0 &=& -3190\:kJmol^{-1} \nonumber -->

<!-- \\ -->

<!-- (\#eq:primary-1) -->

<!-- \\ -->

<!-- \end{eqnarray} -->

<!-- \begin{eqnarray} -->

<!-- (CH_2O)_{106}(NH_3)_{16}(H_3PO_4) &+& 236MnO_2 +472H^+ \rightarrow  \nonumber -->

<!-- \\ -->

<!-- && 106CO_2 + 236Mn^{2+} + 8N_2 + H_3PO_4 + 366H_2O -->

<!-- \\ -->

<!-- \textrm{Free energy, } \Delta G_0 &=& -3050\:kJmol^{-1} \nonumber -->

<!-- \\ -->

<!-- (\#eq:primary-2) -->

<!-- \\ -->

<!-- (CH_2O)_{106}(NH_3)_{16}(H_3PO_4) &+& 94.4HNO_3 \rightarrow  \nonumber -->

<!-- \\ -->

<!-- && 106CO_2 + 55.5N_2 + H_3PO_4 + 177H_2O -->

<!-- \\ -->

<!-- \textrm{Free energy, } \Delta G_0 &=& -3030\:kJmol^{-1} \nonumber -->

<!-- \\ -->

<!-- (\#eq:primary-3) -->

<!-- \\ -->

<!-- (CH_2O)_{106}(NH_3)_{16}(H_3PO_4) &+& 212Fe_2O_3 +848H^+ \rightarrow  \nonumber -->

<!-- \\ -->

<!-- && 106CO_2 + 16NH_3 + H_3PO_4 + 742H_2O + 424Fe^{2+} -->

<!-- \\ -->

<!-- \textrm{Free energy, } \Delta G_0 &=& -1410\:kJmol^{-1} \nonumber -->

<!-- \\ -->

<!-- (\#eq:primary-4) -->

<!-- \\ -->

<!-- (CH_2O)_{106}(NH_3)_{16}(H_3PO_4) &+& 53SO_4^{2-}  \rightarrow  \nonumber -->

<!-- \\ -->

<!-- && 106CO_2 + 16NH_3 + H_3PO_4 + 106H_2O + 53S^{2-} -->

<!-- \\ -->

<!-- \textrm{Free energy, } \Delta G_0 &=& -380\:kJmol^{-1} \nonumber -->

<!-- \\ -->

<!-- (\#eq:primary-5) -->

<!-- \\ -->

<!-- (CH_2O)_{106}(NH_3)_{16}(H_3PO_4) &&  \rightarrow  \nonumber -->

<!-- \\ -->

<!-- && 53CO_2 + 53CH_4 + 16NH_3 + H_3PO_4 -->

<!-- \\ -->

<!-- \textrm{Free energy, } \Delta G_0 &=& -350\:kJmol^{-1} \nonumber -->

<!-- (\#eq:primary-6) -->

<!-- \end{eqnarray} -->

The kinetic rate constant *k~OM~* gives the maximum oxidation rate, which is different for each reactive organic matter type, but the same for each oxidation pathway.

When using `OMModel` 3, the kinetic rate constant is the rate of bacterial growth. The rates are calculated as the product of the rate constant, the concentration and a series of limitation and inhibition factors that scale between 0 and 1.

##### *Limitation and inhibition factors* {.unnumbered}

A common feature of all diagenesis models is that the oxidation rate expression $R_{Ox}$ is a product of up to seven terms:

-   an organic matter reaction rate constant *k~OM~*
-   the organic matter concentration limitation, $F_{OM}$
-   a temperature dependence $F_{Tem}$
-   a microbial biomass factor, $F_{Bio}$
-   a term limitation factor, $F_{TEA}$
-   an inhibition term, $F_{IN}$ and
-   a thermodynamic factor, $F_{T}$
-   a thermodynamic factor, $F_{an}$

<center>

```{=tex}
\begin{equation}
R_{Ox_{i}} = k_{OM}F_{OM}F_{Tem}F_{Bio}F_{TEA{i}}F_{In_{i}}F_{T}F_{an}
\\
(\#eq:factors-1)
\end{equation}
```
</center>

In CANDI-AED, organic matter limitation, $F_{OM}$, is only used with `OMModel` 3. It takes the general form

<center>

```{=tex}
\begin{equation}
F_{OM_{i}} = \frac{OM_{i}} {(OM_{i} + K_{OM})}
\\
(\#eq:factors-2)
\end{equation}
```
</center>

where $OM_i$ is an organic matter species, such as $D_{Hyd}$ or $OAc$ and *K~OM~* is a limiting concentration. This functions like other Monod terms, limiting the reaction rate when the concentration of $OM_{i}$ is low. Using `OMModel` 1 or 2 (the common and simplest options), $F_{OM}$ is effectively 1 and it is assumed that organic matter concentration never limits the reaction rate.

The temperature factor $F_{Tem}$ is rarely employed in diagenesis models, but in a handful of cases it uses a Q~10~ relationship between 2 and 4 (see Fossing et al. 2004 for a clear explanation of how temperature affects reaction rates and Eldridge and Morse (2008) or Reed et al. (2011) for a specific examination of the effect of temperature). The current version of the model has switches built in for the temperature dependence factor, $F_{Tem}$, where values of 1 or 2 turn them off and on. However, implementation and testing of the factors has not been carried out for this version of the model.

The *TEA* factor, $F_{TEA}$, term accounts for the $R_{Ox}$ dependence on the oxidant concentration when the oxidant concentration is low. The $F_{TEA}$ term in Approach 1 is a Monod expression (table \@ref(tab:OMApproach-table), figure \@ref(fig:OMModels-2)), which uses Monod half-saturation constants (*K~Ox~*), and which is chosen because it best reflects laboratory data of bacterially-controlled oxidation reactions (Boudreau and Westrich, 1984, Gaillard and Rabouille 1992). In Approach 2, the $F_{TEA}$ is either 0, 1 or the ratio of *Ox~i~* to *L~Ox~*, depending on the oxidant concentration relative to *L~Ox~*, a specified limiting concentration (table \@ref(tab:OMApproach-table)). We use the notation *K~Ox~* and *L~Ox~* to emphasise that a distinction should be made between the Monod half constants in Approach 1 and the limiting concentrations used in Approach 2; the difference in conceptual representation is not always clear in Approach 2 papers that use the notation *K~Ox~*. We recommend that the user set the parameter `OMModel` to 1, unless they specifically need to emulate an Approach 2 model, such as that by Van Cappellen and Wang (1996).

<center>

```{r OMApproach-table, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
OMApproach_tab <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="OMApproach")
  kable(OMApproach_tab[,1:3],"html", escape = F, align = "c"
        , caption = "Limitation and inhibition equations for the two different approaches (`OMApproach`).",bootstrap_options = "hover")%>%
  
 
   kable_styling(OMApproach_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>% 
    
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "10em" ,color="black",bold = F) %>%
  column_spec(2, width_min = "15em" ,color="black") %>%
  column_spec(3, width_min = "15em" ,color="black") %>%
  # row_spec(1:2, background = 'white')  %>%
 scroll_box(width = "40em", height = "17em",fixed_thead = FALSE)
```

</center>

<br>

<center>

```{r OMModels-2, echo=FALSE, fig.cap="Schematic of examples of how inhibition and limitation scale between 0 and 1 over a range of concentrations. Left: limitation function. Right: inhibtion function.", fig.show="hold", out.width="85%" ,fig.align='center'}
knitr::include_graphics("images/23-sediment_biogeochemistry/OMSetup/FTEAFIn-02.png")
```

</center>

$F_T$ is the thermodynamic factor. The current version of the model includes $F_T$ only for `OMModel` 3, for terminal oxidation reactions and fermentation. This factor uses the ratio of products and reactants, as well as the free energy of each reaction pathway ($\Delta G^0$).

The $F_T$ is a factor that scales between 0, which stops the reaction, and 1, which allows the reaction to proceed, in the same way as $F_{TEA}$ and $F_{In}$. The form of the equation proposed by Jin and Bethke is that shown in equation \@ref(eq:factors-3).The *m* and *χ* are stoichiometric coefficients, *ΔG~ATP~* is the energy required to synthesise one mole of ATP, *R* is the gas constant and *T* is the temperature. In the case where the $F_T$ \< 0 then the factor is set to 0 by convention, otherwise a negative FT would represent the case where microbes carried out a reaction that was unfavourable to their metabolism.

<center>

```{=tex}
\begin{equation}
F_T = 1 - e^{\frac {∆G_r + m ∆G_{ATP}} {χR} } 
(\#eq:factors-3)
\end{equation}
```
</center>

\

*ΔG~r~* is the changing free energy in the system, which controls whether the $F_T$ becomes favourable or not as the environment changes, and is calculated according to the simplified version in equation \@ref(eq:factors-4) (see Jin and Bethke 2005 for a more precise definition):

<center>

```{=tex}
\begin{equation}
    ∆G_r = ∆G^0 + RT ln ( \frac {[redox \ products]^{molar \ ratios} }  { [redox \ reactants]^{molar \ ratios}  } )
\\
(\#eq:factors-4)
\end{equation}
```
</center>

where *ΔG~0~* is the standard state free energy of a redox reaction, calculated as the sum of the energy of the reduction of a *TEA* and the oxidation of an organic substrate. These energies are well constrained by laboratory data, though there may be differences between laboratory and field conditions that create some uncertainty in their use (Bethke et al. 2011). Practically, this equation allows for *ΔG~r~* to be very close to *∆G~0~* when the reactants are at high concentrations relative to the products. Inversely, when the products are at relatively high concentrations, *ΔG~r~* becomes more positive, making $F_T$ more negative. Equation \@ref(eq:factors-3) requires the constants *m* and *χ*, which are parameters that could be measured empirically, however as yet they have not all been measured for every redox pathway for every substrate, and therefore authors have had to estimate some values when using equation \@ref(eq:factors-3) (Dale et al. 2008). It also requires *ΔG~ATP~*, which is around 45 kJ (mol substrate)^-1^ but has varied published values. Around neutral *pH* the conditions for fermentation product consumption by iron reducers, sulphate reducers and methanogens are all thermodynamically favourable, but this changes with *pH* (Bethke et al. 2011). Acidity either promotes or hinders the reaction depending on whether the protons are produced or consumed in the reaction.

$F_{an}$ is a factor that slows down sulfate reduction and methanogenesis. As stated above, the kinetic rate constant *k~OM~* is the same for all redox pathways. The higher-energy pathways are designed to inhibit the lower-energy pathways, but not to occur faster, which is consistent with most other diagenesis models. However, there is an option to slow sulfate reduction and methanogenesis if the user desires. The parameter `f_an` can be set between 0 and 1 to make the $F_{an}$ factor more or less sensitive. Any value greater than 1 makes the factor ineffective (equation \@ref(eq:#factors-4)). The $F_{an}$ factor is calculated as the fraction of $F_{Sul}$ and $F_{Met}$ over the sum of all of the other pathways.

```{=tex}
\begin{equation}
For \ \ R_{O_{2}},\ \ R_{NO_{3}},\ \ R_{MnO_{2}}, \ \ and  \ \ R_{FeOH} \
\\
\\
F_{an (factor)} = 1
\\
\\
For \ \ R_{SO_{4}} \ \ and  \ \ R_{Met} \
\\
\\
F_{an (factor)} =\left\{
  \begin{array}{@{}ll@{}}
    1, & f_{an (parameter)} \: \ge \: 1 \\
    \frac {F_{SO_{4}} + F_{Met}}{f_{an (parameter)} + F_{O_{2}} + F_{NO_{3}} + F_{MnO_{2}} + F_{FeOH} + F_{SO_{4}} + F_{Met}},       & f_{an (parameter)} \lt 1\\
  \end{array}\right.

(#eq:factors-4)

\end{equation}
```
##### *Nitrogen cycling* {#SedNitrogenCycling .unnumbered}

The redox sequence in CANDI-AED is slightly different to the sequence in other diagenesis models: here the nitrogen redox reactions are split into several processes. Most diagenesis models have an overall reaction of $NO_3^-$ oxidising organic matter and producing $N_2$, as shown in figure \@ref(fig:NModel-1). This reaction is less favourable than aerobic oxidation and it is inhibited by $O_2$.

<center>

```{r NModel-1, echo=FALSE, fig.cap="Schematic of the basic nitrogen-organic matter redox process from most other diagenesis models.", fig.show="hold", out.width="35%" ,fig.align='center'}
knitr::include_graphics("images/23-sediment_biogeochemistry/OMSetup/NModel1.png")
```

</center>

In CANDI-AED, there are separate nitrogen reactions for oxidising organic matter:

-   **Denitrousation**: $N_2O$ reduction to $N_2$\
-   **Denitratation**: $NO_3^-$ reduction to $NO_2^-$\
-   **Denitritation**: $NO_2^-$ reduction\
-   **Nitrous denitritation**: $NO_2^-$ reduction to $N_2O$\
-   **DNRA**: (Dissimilatory nitrate reduction to ammonia) $NO_2^-$ reduction to $NH_4^+$

These reactions interact as shown in figure \@ref(fig:NModel-2). In CANDI-AED, the denitrousation reaction is configured to be more favourable than aerobic oxidation, based on the free energy data reported in Lam and Kuypers (2011). Therefore $N_2O$ inhibits aerobic respiration and all lower-energy pathways. Denitratation and all other reactions are less favourable than aerobic oxidation. Denitritation is the term used for the oxidation of $NO_2^-$, the products of which are produced through either nitrous denitritation (to $N_2O$) or DNRA (to $NH_4^+$), depending on the abundance of $NO_2^-$.The rate equations are as shown in table \@ref(tab:NRate-table).

<center>

```{r NModel-2, echo=FALSE, fig.cap="Schematic of the nitrogen-organic matter redox processes used in the CANDI-AED model.", fig.show="hold", out.width="85%" ,fig.align='center'}
knitr::include_graphics("images/23-sediment_biogeochemistry/OMSetup/NModel2.png")
```

<br>

```{r NRate-table, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
NRates_tab <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="NRates")
  kable(NRates_tab[,1:2],"html", escape = F, align = "c"
        , caption = "Rate equations for organic matter oxidation by nitrogen.",bootstrap_options = "hover")%>%
  
 
   kable_styling(NRates_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>% 
    
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "10em" ,color="black",bold = F) %>%
  column_spec(2, width_min = "25em" ,color="black") %>%
  # column_spec(3, width_min = "15em" ,color="black") %>%
  row_spec(1:5, background = 'white')  %>%
 scroll_box(width = "44em", height = "37em",fixed_thead = FALSE)
```

</center>

Along with the nitrogen--organic matter redox reactions, there are sets of nitrogen oxidation reactions, where the oxidant is $O_2$ or $NO_2^-$. These are outlined in more detail in the AED inorganic nitrogen chapter.

##### *Calculation of the primary reactions* {#SedCalculationofthePrimaryReactions .unnumbered}

The chemical reactions in table \@ref(tab:OMChem-table) represent the overall process, however, the model code does not contain them in that form. A description of the numerical solution is presented below in the section [Numerical solution](#SedNumericalSolution) and a summary of the process as it pertains to primary reactions is given here. The model code contains:

**The 'Reaction' function**: the balance of each reaction rate for each species.

**The 'Rates' function**: the reaction rates, as the product of

-   the species,
-   the kinetic constants,
-   the stoichiometric coefficients and
-   the factors.

<center>

```{r OMCalcs-1, echo=FALSE, fig.cap="An illustration of one example of how the reaction rates are calculated in the model code. The balance equation is calculated as the gains and losses from the individual rates, and the rates are calculated from constants or other variables.", fig.show="hold", out.width="100%" ,fig.align='center'}
knitr::include_graphics("images/23-sediment_biogeochemistry/OMSetup/EquationStructure-02.png")
```

</center>

#### Secondary Redox Reactions {#SedSecondaryRedoxReactions}

The reactions of chemical species produced by the primary redox reactions (\@ref(tab:SecondChem-table)) are referred to as secondary redox reactions, and are usually given bimolecular rate laws, which are first order with respect to the oxidant and reductant (\@ref(tab:SecondChem-table)). The rates are controlled by a kinetic constant with units (mmol L)\^-1 y\^-1.

<center>

```{r SecondChem-table, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
SecondChem_tab <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="SecondChem (3)")
  kable(SecondChem_tab[,1:5],"html", escape = F, align = "c"
        , caption = "Chemical reactions for secondary redox reactions.",bootstrap_options = "hover")%>%
   kable_styling(SecondChem_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 11) %>% 
    
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "15em" ,color="black",bold = T) %>%
  column_spec(2, width_min = "3em" ,color="black") %>%
  column_spec(3, width_min = "12em" ,color="black") %>%
  column_spec(4, width_min = "3em" ,color="black") %>%
  column_spec(5, width_min = "12em" ,color="black") %>%
  row_spec(1:5, background = 'white')  %>%
 scroll_box(width = "770px", height = "45em",fixed_thead = FALSE)
```

</center>

##### *Calculation of the secondary reactions* {#SedCalculationoftheSecondaryReactions .unnumbered}

As with the [primary reactions](SedCalculationofthePrimaryReactions), the secondary redox reactions are present in the code as the combination of

**The 'Reaction' function**: the sum of the rates of change of each species.

**The 'Rates' function**: the reaction rates, as the product of

-   the species,
-   the kinetic constants,
-   the stoichiometric coefficients and
-   the factors.

<center>

```{r OMCalcs-2, echo=FALSE, fig.cap="An illustration of one example of how the reaction rates are calculated in the model code. The balance equation is calculated as the gains and losses from the individual rates, and the rates are calculated from constants or other variables.", fig.show="hold", out.width="100%" ,fig.align='center'}
knitr::include_graphics("images/23-sediment_biogeochemistry/OMSetup/EquationStructure-03.png")
```

</center>

#### Equilibrium Geochemistry {#EquilibriumGeochemistry}

##### *pH* {.unnumbered}

The *pH* is initialised at a fixed value of around 8, which the user does not control. Setting the parameter `rxn_mode` to 0 or 4 keeps the *pH* constant throughout the simulation. Setting `rxn_mode` to 1, 2 or 3 updates the *pH* calculation at each time step. The *pH* is calculated as the sum of all charged species, where any unbalanced charge is calculated as 'charge balance' ($ubalchg$), which is assumed to be H^+^. The charge balance is at each time step is solved as a state variable, which is subject to advection, diffusion and bioturbation reactions.

The *pH* is used in the calculation of $PO_4^{3-}$ adsorption if the parameter `ads_use_pH` is set to `TRUE`.

*pH* is not used in any of the other calculations.

##### *Mineral precipitation and ageing* {#SedMineralPrecipitationAndAgeing .unnumbered}

In CANDI-AED, dissolved species can combine to form precipitated minerals: $MnO_{2A}$, $Fe(OH)_{3A}$, $FeS$, $MnCO_3$, $FeCO_3$, and $CaCO_3$. Another process is the ageing of amorphous $MnO_{2A}$ and $Fe(OH)_{3A}$ to the unreactive, crystalline phases $MnO_{2B}$ and $Fe(OH)_{3B}$. The third type of process is ageing process of $FeS$ to $FeS_2$ (adding a further S atom to the compound). The reactions rates are outlined in tables \@ref(tab:Geochem-table1) to \@ref(tab:Geochem-table3).

The form of the rate is dependent on the complexity of the reaction mode chosen, using the parameter `rxn_mode`. The simplest option is for the reactions not to occur. The next option is for the reaction to proceed as a biomolecular rate law, in a similar way to the secondary redox reactions, by selecting `rxn_mode` = 2.

There are also two methods that set the precipitation rate based on the relative concentrations of dissolved and precipitated species. One method uses the ion activity product (*IAP*), which is a state variable and is calculated at each depth and timestep. The calculation method is based on that from Tufano et al. (2009). This method is generally selected with `rxn_mode` = 3 (see tables \@ref(tab:Geochem-table1) to \@ref(tab:Geochem-table3)). The other method is the "omega" method, using the equations from Van Cappellen and Wang (1996), generally by using `rxn_mode` = 4 (see tables \@ref(tab:Geochem-table1) to \@ref(tab:Geochem-table3)).

<center>

```{r Geochem-table1, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
MnOFeO_tab    <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="MnOFeO")
MnOFeO_Groups <- unique(MnOFeO_tab$Group)
  kable(MnOFeO_tab[,1:4],"html", escape = F, align = "c"
        , caption = "Reactions and rate equations for precipitation and ageing for MnO~2~ and Fe(OH)~3~",bootstrap_options = "hover")%>%
    
  pack_rows(MnOFeO_Groups[1],
            min(which(MnOFeO_tab$Group == MnOFeO_Groups[1])),
            max(which(MnOFeO_tab$Group == MnOFeO_Groups[1])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(MnOFeO_Groups[2],
            min(which(MnOFeO_tab$Group == MnOFeO_Groups[2])),
            max(which(MnOFeO_tab$Group == MnOFeO_Groups[2])),
            background = '#ebebeb'
            ,color="black") %>%
    
   kable_styling(MnOFeO_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>% 
    
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "10em" ,color="black",bold = F) %>%
  column_spec(2, width_min = "15em" ,color="black") %>%
  column_spec(3, width_min = "3em"  ,color="black") %>%
  column_spec(4, width_min = "15em" ,color="black") %>%
  # row_spec(1:4, background = 'white')  %>%
 scroll_box(width = "48em", height = "30em",fixed_thead = FALSE)
```

</center>

$MnO_{2A}$ ageing

```{=tex}
\begin{equation}
MnO_{2A} \rightarrow MnO_{2B}
\end{equation}
```
::: {.panelset .sideways}
::: panel
[rxn_mode = 0]{.panel-name}

```{=tex}
\begin{equation}
 R_{MnO_{2B}ppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 1]{.panel-name}

```{=tex}
\begin{equation}
 R_{MnO_{2B}ppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 2]{.panel-name}

```{=tex}
\begin{equation}
 R_{MnO_{2B}ppt} = k_{MnO_{2B}ppt}  Mn_{O_{2A}}
\end{equation}
```
:::

::: panel
[rxn_mode = 3]{.panel-name}

```{=tex}
\begin{equation}
 R_{MnO_{2B}ppt} = k_{MnO_{2B}ppt}  Mn_{O_{2A}}
\end{equation}
```
:::

::: panel
[rxn_mode = 4]{.panel-name}

```{=tex}
\begin{equation}
 R_{MnO_{2B}ppt} = k_{MnO_{2B}ppt}  Mn_{O_{2A}}
\end{equation}
```
:::
:::

\
\
<!-- _______________________________________________________________________________ -->

$MnO_{2A}$ precipitation

```{=tex}
\begin{equation}
2Mn^{2+} + O_2 + 4HCO_3^- \rightarrow  2MnO_{2A} + 4CO_2 +2H_2O
\end{equation}
```
::: {.panelset .sideways}
::: panel
[rxn_mode = 0]{.panel-name}

```{=tex}
\begin{equation}
 R_{MnO_{2A}ppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 1]{.panel-name}

```{=tex}
\begin{equation}
 R_{MnO_{2A}ppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 2]{.panel-name} Continuous

```{=tex}
\begin{equation}
 R_{MnO_{2A}ppt} = k_{MnO_{2A}ppt}  Mn^{2+}
\end{equation}
```
:::

::: panel
[rxn_mode = 3]{.panel-name} IAP method

```{=tex}
\begin{equation}

 R_{MnO_{2A}ppt} = \left\{
  \begin{array}{@{}ll@{}}
    k_{MnO_{2A}ppt}  MnO_{2A}, & IAP_{MnO2} = 0 \\
    k_{MnO_{2A}ppt} \left( 1 - \frac {1}{IAP_{MnO2}} \right) , & IAP_{MnO2} \lt 1 \\
    k_{MnO_{2A}ppt} \left( 1 - \frac {IAP_{MnO2}}{1} \right)   & IAP_{MnO2} \ge 1
  \end{array}\right.
\end{equation}
```
:::

::: panel
[rxn_mode = 4]{.panel-name} Continuous

```{=tex}
\begin{equation}
 R_{MnO_{2A}ppt} = k_{MnO_{2A}ppt}  Mn^{2+}
\end{equation}
```
:::
:::

<!-- _______________________________________________________________________________ -->

\
\

$Fe(OH)_{3A}$ ageing

```{=tex}
\begin{equation}
Fe(OH)_{3A} \rightarrow Fe(OH)_{3B}
\end{equation}
```
::: {.panelset .sideways}
::: panel
[rxn_mode = 0]{.panel-name}

```{=tex}
\begin{equation}
 R_{Fe(OH)_{3A}ppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 1]{.panel-name}

```{=tex}
\begin{equation}
R_{Fe(OH)_{3A}ppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 2]{.panel-name}

```{=tex}
\begin{equation}
 R_{Fe(OH)_{3A}ppt} = k_{Fe(OH)_{3A}ppt}  Fe(OH)_{3A}
\end{equation}
```
:::

::: panel
[rxn_mode = 3]{.panel-name}

```{=tex}
\begin{equation}
 R_{Fe(OH)_{3A}ppt} = k_{Fe(OH)_{3A}ppt} Fe(OH)_{3A}
\end{equation}
```
:::

::: panel
[rxn_mode = 4]{.panel-name}

```{=tex}
\begin{equation}
 R_{Fe(OH)_{3A}ppt} = k_{Fe(OH)_{3A}ppt}  Fe(OH)_{3A}
\end{equation}
```
:::
:::

\
\
\
\
<!-- _______________________________________________________________________________ --> $Fe(OH)_{3A}$ precipitation

```{=tex}
\begin{equation}
4Fe^{2+} + O_2 +8HCO_3^- + 2H_2O \rightarrow 4Fe(OH)_{3A} + 8CO_2 
\end{equation}
```
::: {.panelset .sideways}
::: panel
[rxn_mode = 0]{.panel-name}

```{=tex}
\begin{equation}
 R_{Fe(OH)_{3A}ppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 1]{.panel-name}

```{=tex}
\begin{equation}
 R_{Fe(OH)_{3A}ppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 2]{.panel-name} Continuous

```{=tex}
\begin{equation}
 R_{Fe(OH)_{3A}ppt} = k_{Fe(OH)_{3A}ppt}  Fe^{2+}
\end{equation}
```
:::

::: panel
[rxn_mode = 3]{.panel-name} IAP method

```{=tex}
\begin{equation}

 R_{Fe(OH)_{3A}ppt} = \left\{
  \begin{array}{@{}ll@{}}
    k_{Fe(OH)_{3A}ppt}  Fe(OH)_{3A}, & IAP_{Fe(OH)_3} = 0 \\
    k_{Fe(OH)_{3A}ppt} \left( 1 - \frac {1}{IAP_{Fe(OH)_3}} \right) , & IAP_{Fe(OH)_3} \lt 1 \\
    k_{Fe(OH)_{3A}ppt} \left( 1 - \frac {IAP_{Fe(OH)_3}}{1} \right)   & IAP_{Fe(OH)_3} \ge 1
  \end{array}\right.
\end{equation}
```
:::

::: panel
[rxn_mode = 4]{.panel-name} Continuous

```{=tex}
\begin{equation}
 R_{Fe(OH)_{3A}ppt} = k_{Fe(OH)_{3A}ppt}  Fe^{2+}
\end{equation}
```
:::
:::

\
\
\
\
\
\
\
\
\
\
\
\
\
\
<!-- _______________________________________________________________________________ -->\
\
\
\
\
\
\
<!-- _______________________________________________________________________________ -->

<br>

<center>

```{r Geochem-table2, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
FeS_tab <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="FeS")
# View(sediment_VarsTable)
  kable(FeS_tab[,1:4],"html", escape = F, align = "c"
        , caption = "Reactions and rate equations for precipitation and ageing for FeS and FeS~2~",bootstrap_options = "hover")%>%
   kable_styling(FeS_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>% 
    
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "10em" ,color="black",bold = F) %>%
  column_spec(2, width_min = "15em" ,color="black") %>%
  column_spec(3, width_min = "3em" ,color="black") %>%
  column_spec(4, width_min = "15em" ,color="black") %>%
  row_spec(1:4, background = 'white')  %>%
 scroll_box(width = "48em", height = "30em",fixed_thead = FALSE)
```

</center>

<!-- _______________________________________________________________________________ -->

$FeS$ transformation to $FeS_2$

```{=tex}
\begin{equation}
FeS + H_2S \rightarrow FeS_2 + H_2 \\
R_{Pyrite} = k_{Pyrite} [FeS] [H_2S]
\end{equation}
```
<!-- _______________________________________________________________________________ -->

$FeS$ precipitation

```{=tex}
\begin{equation}
Fe^{2+} + H_2S  \rightarrow FeS + 2H^+
\end{equation}
```
::: {.panelset .sideways}
::: panel
[rxn_mode = 0]{.panel-name}

```{=tex}
\begin{equation}
 R_{FeS} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 1]{.panel-name}

```{=tex}
\begin{equation}
 R_{FeS} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 2]{.panel-name} Continuous

```{=tex}
\begin{equation}
 R_{FeS} = k_{FeSppt}  [Fe^{2+}][HS^-]
\end{equation}
```
:::

::: panel
[rxn_mode = 3]{.panel-name} IAP method

```{=tex}
\begin{equation}
  R_{FeSppt} = k_{FeSppt} \left( 1- \frac{K_{spFeS}}{IAP_{FeS}} \right) \\
  R_{FeSdiss} = - k_{FeSppt} \left( 1- \frac{IAP_{FeS}}{K_{spFeS}} \right) 
  
\end{equation}
```
:::

::: panel
[rxn_mode = 4]{.panel-name} Omega method

```{=tex}
\begin{equation}
R_{FeSppt} = k_{FeSppt} \delta_{FeS} (\Omega_{FeS} - 1)
\\
R_{FeSdiss} = k_{FeSdiss} \delta_{-FeS} (1 - \Omega_{FeS} )
\\
\Omega_{FeS} = \left( \frac {[Fe^{2+}]{HS^-} }    {[H^+] K_{FeS}} \right)
\\
If \Omega_{FeS} \gt 1 \
\delta_{FeS} = 1; \delta_{-FeS} = 0\
\\
If \Omega_{FeS} \le 1 \
\delta_{FeS} = 0; \delta_{-FeS} = 1
\end{equation}
```
:::
:::

<!-- _______________________________________________________________________________ -->

\
\
\
\
\
\

<!-- _______________________________________________________________________________ -->

$MnCO_3$ precipitation

```{=tex}
\begin{equation}
Mn^{2+} + CO_3^{2-} \rightarrow  MnCO_3
\end{equation}
```
::: {.panelset .sideways}
::: panel
[rxn_mode = 0]{.panel-name}

```{=tex}
\begin{equation}
 R_{Rodppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 1]{.panel-name}

```{=tex}
\begin{equation}
 R_{Rodppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 2]{.panel-name} Continuous

```{=tex}
\begin{equation}
 R_{Rodppt} = k_{Rodppt}  [Mn^{2+}][CO_3^{2-}]
\end{equation}
```
:::

::: panel
[rxn_mode = 3]{.panel-name} IAP method

```{=tex}
\begin{equation}
  R_{Rodppt} = k_{Rodppt} \left( 1- \frac{K_{spRod}}{IAP_{Rod}} \right) \\
  R_{Rodppt} = - k_{Rodppt} \left( 1- \frac{IAP_{Rod}}{K_{spRod}} \right) 
  
\end{equation}
```
:::

::: panel
[rxn_mode = 4]{.panel-name} Omega method

```{=tex}
\begin{equation}
R_{Rodppt} = k_{Rodppt} \delta_{Rod} (\Omega_{Rod} - 1)
\\
R_{Roddiss} = k_{Roddiss} \delta_{-Rod} (1 - \Omega_{Rod} )
\\
\Omega_{Rod} = \left( \frac {[Mn^{2+}][CO_3^{2-}] }    {K_{Rod}} \right)
\\
If \  \  \Omega_{Rod} \gt 1 \
\delta_{Rod} = 1; \delta_{-Rod} = 0\
\\
If \  \  \Omega_{Rod} \le 1 \
\delta_{Rod} = 0; \delta_{-Rod} = 1
\end{equation}
```
:::
:::

<!-- _______________________________________________________________________________ -->

\
\
\
\
\
\

<!-- _______________________________________________________________________________ -->

$FeCO_3$ precipitation

```{=tex}
\begin{equation}
Fe^{2+} + CO_3^{2-} \rightarrow  FeCO_3
\end{equation}
```
::: {.panelset .sideways}
::: panel
[rxn_mode = 0]{.panel-name}

```{=tex}
\begin{equation}
 R_{Sidppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 1]{.panel-name}

```{=tex}
\begin{equation}
 R_{Sidppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 2]{.panel-name} Continuous

```{=tex}
\begin{equation}
 R_{Sidppt} = k_{Sidppt}  [Fe^{2+}][CO_3^{2-}]
\end{equation}
```
:::

::: panel
[rxn_mode = 3]{.panel-name} IAP method

```{=tex}
\begin{equation}
  R_{Sidppt} = k_{Sidppt} \left( 1- \frac{K_{spSid}}{IAP_{Sid}} \right) \\
  R_{Sidppt} = - k_{Sidppt} \left( 1- \frac{IAP_{Sid}}{K_{spSid}} \right) 
  
\end{equation}
```
:::

::: panel
[rxn_mode = 4]{.panel-name} Omega method

```{=tex}
\begin{equation}
R_{Sidppt} = k_{Sidppt} \delta_{Sid} (\Omega_{Sid} - 1)
\\
R_{Siddiss} = k_{Siddiss} \delta_{-Sid} (1 - \Omega_{Sid} )
\\
\Omega_{Sid} = \left( \frac {[Fe^{2+}][CO_3^{2-}] }    {K_{Sid}} \right)
\\
If \  \  \Omega_{Sid} \gt 1 \
\delta_{Sid} = 1; \delta_{-Sid} = 0\
\\
If \  \  \Omega_{Sid} \le 1 \
\delta_{Sid} = 0; \delta_{-Sid} = 1
\end{equation}
```
:::
:::

<!-- _______________________________________________________________________________ -->

\
\
\
\
\
\

<!-- _______________________________________________________________________________ -->

$CaCO_3$ precipitation

```{=tex}
\begin{equation}
Ca^{2+} + CO_3^{2-} \rightarrow  CaCO_3
\end{equation}
```
::: {.panelset .sideways}
::: panel
[rxn_mode = 0]{.panel-name}

```{=tex}
\begin{equation}
 R_{Calppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 1]{.panel-name}

```{=tex}
\begin{equation}
 R_{Calppt} = 0
\end{equation}
```
:::

::: panel
[rxn_mode = 2]{.panel-name} Continuous

```{=tex}
\begin{equation}
 R_{Calppt} = k_{Calppt}  [Ca^{2+}][CO_3^{2-}]
\end{equation}
```
:::

::: panel
[rxn_mode = 3]{.panel-name} IAP method

```{=tex}
\begin{equation}
  R_{Calppt} = k_{Calppt} \left( 1- \frac{K_{spCal}}{IAP_{Cal}} \right) \\
  R_{Calppt} = - k_{Calppt} \left( 1- \frac{IAP_{Cal}}{K_{spCal}} \right) 
  
\end{equation}
```
:::

::: panel
[rxn_mode = 4]{.panel-name} Omega method

```{=tex}
\begin{equation}
R_{Calppt} = k_{Calppt} \delta_{Cal} (\Omega_{Cal} - 1)
\\
R_{Caldiss} = k_{Caldiss} \delta_{-Cal} (1 - \Omega_{Cal} )
\\
\Omega_{Cal} = \left( \frac {[Ca^{2+}][CO_3^{2-}] }    {K_{Cal}} \right)
\\
If \  \ \Omega_{Cal} \gt 1 \
\delta_{Cal} = 1; \delta_{-Cal} = 0\
\\
If \  \  \Omega_{Cal} \le 1 \
\delta_{Cal} = 0; \delta_{-Cal} = 1
\end{equation}
```
:::
:::

<!-- _______________________________________________________________________________ -->

<br>

<center>

```{r Geochem-table3, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
CO3_tab <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="CO3")
CO3_Groups <- unique(CO3_tab$Group)
  kable(CO3_tab[,1:4],"html", escape = F, align = "c"
        , caption = "Reactions and rate equations for precipitation and ageing for carbonate compounds.",bootstrap_options = "hover")%>%
   pack_rows(CO3_Groups[1],
            min(which(CO3_tab$Group == CO3_Groups[1])),
            max(which(CO3_tab$Group == CO3_Groups[1])),
            background = '#ebebeb'
            ,color="black") %>%
   pack_rows(CO3_Groups[2],
            min(which(CO3_tab$Group == CO3_Groups[2])),
            max(which(CO3_tab$Group == CO3_Groups[2])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(CO3_Groups[3],
            min(which(CO3_tab$Group == CO3_Groups[3])),
            max(which(CO3_tab$Group == CO3_Groups[3])),
            background = '#ebebeb'
            ,color="black") %>%
    
    
   kable_styling(CO3_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>% 
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "10em" ,color="black",bold = F) %>%
  column_spec(2, width_min = "15em" ,color="black") %>%
  column_spec(3, width_min = "3em" ,color="black") %>%
  column_spec(4, width_min = "15em" ,color="black") %>%
  row_spec(1:4, background = 'white')  %>%
 scroll_box(width = "48em", height = "30em",fixed_thead = FALSE)
```

</center>

<br>

The parameter `rxn_mode` sets different configurations of precipitation and *pH* calculations. A summary of the options is presented in \@ref(tab:Geochem-table4).

<center>

```{r Geochem-table4, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
rxn_mode_tab <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="rxn_mode")
  kable(rxn_mode_tab[,1:9],"html", escape = F, align = "c"
        , caption = "Geochemical configurations from the `rxn_mode` parameter",bootstrap_options = "hover")%>%
   kable_styling(rxn_mode_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>% 
    
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "12em" ,color="black",bold = T) %>%
  column_spec(2, width_min = "12em" ,color="black") %>%
  column_spec(3, width_min = "12em" ,color="black") %>%
  column_spec(4, width_min = "12em" ,color="black") %>%
  column_spec(9, width_min = "12em" ,color="red", italic = T) %>%
  row_spec(1:4, background = 'white')  %>%
 scroll_box(width = "48em", height = "25em",fixed_thead = FALSE)
```

</center>

##### *Adsorption* {.unnumbered}

IN CANDI-AED, $NH_4^+$ adsorption is simulated as precipitating onto solid particles, including organic matter and clay, following the method of Van Cappellen and Wang (1996). It is assumed that there are unlimited adsorption sites for $NH_4^+$, and the relative amounts of $NH_4^+$ were assigned by the adsorption constant `knh4p`, which was initially based on Mackin and Aller (1984). A simplified version of the adsorption equations adsorption can be found in the equation set below (equation \@ref(eq:ppt-1)). The full equations correct for porosity and solid volume when transferring between particulate and dissolved phases. $NH_4^+$ adsorption defaults ot being switched on, but may be switched off by setting the parameter `NH4AdsorptionModel` to 0.

<center>

```{=tex}
\begin{eqnarray}

NH_{4 _{Total}}^+       = NH_{4}^+ + NH_{4_{Adsorbed}}
\\
\\
NH_{4 _{Dissolved}}^+   = \frac {1}{K_{NH_4p}}
\\
\\
NH_{4 _{Particulate}}^+ = 1 - NH_{4 _{Dissolved}}^+ F_{Sal}
\\
\\
NH_{4 _{Particulate}}^+ = NH_{4_{Total}}^+ NH_{4 _{Particulate}}^+
\\
\\
NH_{4 _{Dissolved}}^+   = NH_{4_{Total}}^+ - NH_{4 _{Particulate}}^+
(\#eq:ppt-1)
\end{eqnarray}
```
</center>

A similar approach is taken for dissolved organic carbon, nitrogen and phosphrous. The $DON_R$ and $DON_S$ are adsorbed to the sediment with the adsorption constant *K~DOMP~*, and $DOP_R$ and $DOP_S$ adsorbed using the same adsorption constant *K~DOMP~*. The value for *K~DOMP~* defaults to 1.4 and the model defaults to being switched on (`DOMAdsorptionModel` \> 0).

The model uses the parameter `PO4AdsorptionModel` to switch the type of adsorption equations from the AED $PO_4^{3-}$ adsorption library. If it is not set, the model defaults to method 1. The equations for `PO4AdsorptionModel` 1 are given in equation set \@ref(eq:ppt-2). $PO_4^{3-}$ is adsorbed onto the surfaces of iron oxide particles, both the reactive $Fe(OH)_{3A}$ and unreactive $Fe(OH)_{3B}$ fractions. During a simulation, low oxygen and nitrate concentrations in the sediment cause iron oxide reduction, which dissolves the solid particles and releases the dissolved $PO_4^{3-}$ to the porewater.

If `PO4AdsorptionModel` 2 is selected, a set of pH-dependent reactions are used. A further parameter, `ads_use_pH` is required to be set to `TRUE` if the user wants the model to respond to *pH* dynamically, or to `FALSE` if it should be less sensitive to *pH*. If `PO4AdsorptionModel` 3 is selected, a set of equations is used that rely on salinity and temperature. It is suggested to use `PO4AdsorptionModel` = 1 for most sediment model applications.

```{=tex}
\begin{eqnarray}

PO_{4_{Particulate}}^{3-} = \frac {K_{PO_4p} (Fe(OH)_{3A} + Fe(OH)_{3B})}
{1 + K_{PO_4p} (Fe(OH)_{3A} + Fe(OH)_{3B}) } (PO_{4 _{Total}}^{3-})
\\
\\
PO_{4_{Dissolved}}^{3-}   = \frac {1} {1 + K_{PO_4p} (Fe(OH)_{3A} + Fe(OH)_{3B}) } (PO_{4 _{Total}}^{3-})
(\#eq:ppt-2)
\end{eqnarray}
```
#### Biological dynamics

Animals and plants are included in determining the transport and biogeochemistry of the sediment.

##### *Benthic macrofauna* {.unnumbered}

###### *Bioturbation* {.unnumbered}

Bioturbation is a mixing process caused by benthic infauna in the upper layers of the sediment that affects both solids and porewater (Berg et al. 2001). The model uses a diffusion coefficient (bio-diffusivity, $D_B$) that varies with depth ($D_B[z_s]$) as a two layer function or a Gaussian decrease (Boudreau 1996). Bioirrigation is a non-local mixing process in the upper layers of the sediment, caused by animal burrows leading to enhanced diffusion of porewater (figure \@ref(fig:Phys3)).

```{r Phys3, echo = FALSE, out.width='60%', class = "text-image",fig.show='hold',fig.align='center', fig.cap = "Bioturbation is the mixing of solids and porewater in the upper layers by benthic animals. Irrigation causes extra mixing of porewater through their burrows."}

knitr::include_graphics("images/23-sediment_biogeochemistry/AdvectionDispersion/Biot-07.png")

```

Bioturbation is computed to vary with depth below the sediment-water interface, based on selection of an appropriate model, $\Theta_{\text{imix}}^{sdg}$, via the switch `imix` (equation \@ref(eq:biodynamics-1), figure \@ref(fig:Init-3)), such that:

```{=tex}
\begin{equation}
  \Theta_{\text{imix}}^{sdg}=\left\{
  \begin{array}{@{}ll@{}}
    0, & \text{exponential decay curve}\\
    1, & \text{piece-wise two-layer option}\\
    2  & \text{exponential decay curve}\\
  \end{array}\right.
\end{equation}
(\#eq:biodynamics-1)
```
<br>

```{r Init-3, echo=FALSE, fig.cap="Depth profile representing two possible bioturbation options and the functions of the parameters `DB0`, `x1` and `x2`. ",fig.show='hold',fig.align='center' ,out.width = '80%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Initialisation/Bioturb-02.png")
```

The rate of bioturbation is set with the parameter `DB0`, whereby the bio-diffusivity is computed as $D_{B}=D_{B}[D_{B_0},z_s]$. If the user selects `imix` switch 0 or 2, then they use the parameter `xs` to set the shape of the exponential decay curve. If the user selects `imix` switch 1, then they use parameters `x1` and `x2` to set the two depths between which bioturbation goes at the full rate and drops to zero:

::: {.panelset .sideways}
::: panel
[imix = 0]{.panel-name} Bio-diffusivity:

```{=tex}
\begin{equation}
 D_B[z_s] = D_{B_0} e^{-\frac {z_s^2}{2 xs^2}}  
\end{equation}
```
:::

::: panel
[imix = 1]{.panel-name} Bio-diffusivity:

```{=tex}
\begin{equation}
  D_B[z_s] =\left\{
  \begin{array}{@{}ll@{}}
    D_{B_0}, & z_s \: \le \: x_1 \\
    D_{B_0} \: (\frac  {x_2 - z_s} {x_2 - x_1}),       & z_s \gt x_1 \text{ and } z_s \le x_2\\
    0        & z_s \gt x_2 \\
  \end{array}\right.
\end{equation}
```
:::

::: panel
[imix = 2]{.panel-name} Bio-diffusivity:

$D_B[z_s] = D_{B_0} \: e^{(- \frac{z_s}{xs} )}$
:::
:::

A bioturbation profile is output in the file `Depths.sed`, so the user can examine the shape of the profile based on the parameter settings chosen.

###### *Bioirrigation* {.unnumbered}

Bioirrigation has only one option, whereby the parameter `alpha0`, $\alpha_0$, sets the rate of irrigation and `xirrig`, $x_{irrig}$, sets the depth at which irrigation decays (equation \@ref(eq:biodynamics-1), figure \@ref(fig:Init-4)), according to:

::: {.panelset .sideways}
::: panel
[alpha]{.panel-name} Bio-irrigation:

```{=tex}
\begin{equation}
  \alpha[z_s] = \left\{
  \begin{array}{@{}ll@{}}
     \alpha_0 , & z_s \: \le \: x_{irrig} \\
     \alpha_0 \: e^{-1.5 \: (z_s - x_{irrig})},       & z_s \gt x_{irrig} \\
  \end{array}\right.
\end{equation}
(\#eq:biodynamics-2)
```
:::
:::

<br>

```{r Init-4, echo=FALSE, fig.cap="Example of a depth profile that sets the rate of irrigation, and the parameters `alpha_0` and `x_irrig`.",fig.show='hold',fig.align='center', out.width = '50%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Initialisation/Irrig-01.png")
```

###### *Salinity inhibition* {.unnumbered}

A salinity factor, $F_{Sal}$, inhibits processes at hypersaline concentrations. $F_{Sal}$ ranges from 1 to 0, as set between salinity concentration parameters `Sal1` and `Sal2` (PSU), according to equation \@ref(eq:sal-1).

```{=tex}
\begin{equation}
  F_{Sal} =\left\{
  \begin{array}{@{}ll@{}}
    1, & Salinity \: \le \: Sal_1 \\
    \frac {Sal_2 - Salinity} {Sal_2 - Sal_1},       & Salinity \gt Sal_1 \\
    0        & Salinity \gt Sal_2 \\
  \end{array}\right.
(\#eq:sal-1)
\end{equation}
```
The effect of the scaling is shown in figure \@ref(fig:PhysChem1).

```{r PhysChem1, echo = FALSE, out.width='50%', class = "text-image",fig.show='hold',fig.align='center', fig.cap = "Diagram of the salinity factor $F_{Sal}$, scaling between 1 and 0 between the parameters. In this diagram `Sal~1~` is 40 and `Sal~2~` is 70."}
knitr::include_graphics("images/23-sediment_biogeochemistry/PhysChemInteractions/FSal-02.png")
```

$F_{Sal}$ is applied to biota-driven mixing (see below) and to the nitrogen redox and adsorption reactions. It is multiplied by the following reactions:

-   *R~N2O~*
-   *R~NO2~*
-   *R~NO3~*
-   *R~NH4ox~*
-   *R~NH4NO2~*
-   *R~NO2O2~*
-   *NH~4~^+^ adsorption*

The term affects all rates equally and so does not favour the concentration of any of the species. However, it slows down the transformation of each species into the other, and increased the effect of transport reactions and depth-driven concentration gradients.

###### *H~2~S inhibition* {.unnumbered}

High $H_{2}S$ concentration is set to inhibit biota-driven mixing and nitrogen redox reactions, on the basis that $H_{2}S$ is toxic to organisms. An $H_{2}S$ factor, $F_{Sul}$, is an inhibition factor similar to other inhibition factors in the sediment model.

```{=tex}
\begin{eqnarray}

F_{Sul} = \frac {K_{H_{2}S}} {(K_{H_{2}S} + H_{2}S)}
(\#eq:sdg-6)
\end{eqnarray}
```
As with $F_{Sal}$, $F_{Sul}$ scales between 1 and 0, however, as a decay rather than a stepped function (Figure \@ref(fig:PhysChem2)). The value for `KH2S` can be found in the physical parameters table below In this simulation, the same value for `KH2S` was used in all zones.

```{r PhysChem2, echo = FALSE, out.width='50%', class = "text-image",fig.show='hold',fig.align='center', fig.cap = "Diagram of the $H_2S$ factor, $F_{Sul}$, scaling from 1 to 0 as $H_2S$ concentration decreases. "}

knitr::include_graphics("images/23-sediment_biogeochemistry/PhysChemInteractions/FSul-02.png")

```

\

##### *Benthic flora* {.unnumbered}

CANDI-AED has two types of benthic plants that can affect $O_2$ concentrations and physical mixing. They have optional links to the water column models.

###### *Microphytobenthos* {.unnumbered}

Microphytobenthos ($MPB$) are small units of phytoplankton that grow on the sediment surface. They are subject to physical transport processes, in the same was as other sediment solids. $MPB$ concentration increases by growing (rather than reacting), which is the difference between its rates of growth (gross primary production, *R~gpp~*) and respiration (*R~rsp~*) (equation \@ref(eq:MAG-1)).

```{=tex}
\begin{equation}
\frac {\delta MPB} {\delta T} = R_{gpp} - R_{rsp}
(\#eq:MAG-1)
\end{equation}
```
*R~gpp~* and *R~rsp~* are calculated according to equations \@ref(eq:MAG-2) and \@ref(eq:MAG-3). `MPBdepth` and `fgpp_sflux` are CANDI-AED parameters, set for each zone in the parameters input file. *MPBG* and *MPBR* are growth and respiration rate variables linked to the water column model. If they are not present in the water column model, *MPBG* and *MPBR* default to zero, which in turn sets *R~gpp~* and *R~rsp~* to zero. From the sediment-water interface down to `MPBdepth`, the $MPB$ both grow and respire, and below `MPBdepth`, any $MPB$ that are mixed down to those layers cease to grow, due to a lack of light, but continue to respire at a reduced rate.

```{=tex}
\begin{equation}
  R_{gpp} =\left\{
  \begin{array}{@{}ll@{}}
    0, & depth \: \gt \: MPBdepth \\
    \frac{MPBG \times \delta depth} {depth} \times (1 - fgpp - sflux) & depth \le MPBdepth \\
  \end{array}\right.
(\#eq:MAG-2)
\end{equation}
```
```{=tex}
\begin{equation}
  R_{rsp} =\left\{
  \begin{array}{@{}ll@{}}
    \frac {365}{7} \times MPB, & depth \: \gt \: MPBdepth \\
    \frac{MPBR \times \delta depth} {depth} \times (1 - fgpp - sflux) & depth \le MPBdepth \\
  \end{array}\right.
(\#eq:MAG-3)
\end{equation}
```
*R~gpp~* is a source of $O_2$ in the $O_2$ balance equation. *R~rsp~* consumes $O_2$ and produces organic matter.

```{r PhysChem3, echo = FALSE, out.width='50%', class = "text-image",fig.show='hold',fig.align='center', fig.cap = "Schematic of the parameterisation of $O_2$ injection into the sediment where $MPB$ are respiring."}
knitr::include_graphics("images/23-sediment_biogeochemistry/PhysChemInteractions/O2Inject-11.png")

```

###### *Aquatic macrophytes* {.unnumbered}

CANDI-AED also includes aquatic macrophytes, such as seagrass or filamentous algae. In a similar way to $MPB$, the roots of macrophtes inject $O_2$ into the upper layers of the sediment. The rate of $O_2$ production, *R~rootsO2~*, is calculated according to \@ref(eq:MAG-6).

```{=tex}
\begin{equation}
  R_{RootsO_{2}} =\left\{
  \begin{array}{@{}ll@{}}
    0, & depth \: \gt \: RTDP \\
    \frac {rto2 \times \delta depth}{depth}   & depth \le RTDP \\
  \end{array}\right.
(\#eq:MAG-6)
\end{equation}
```
The term `RTDP` is a CANDI-AED parameter, set for each zone. The term *rto2* is the $O_2$ excretion rate parameter. It may be sourced from the water column model if the models are linked, or it may be calculated within CANDI-AED if the models are not linked. The options for having the root $O_2$ settins linked or not linked are laid out in table \@ref(tab:MAG-table).

<center>

```{r MAG-table, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
MAG_tab <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="MAG")
MAG_Groups <- unique(MAG_tab$Group)
  kable(MAG_tab[,1:3],"html", escape = F, align = "c"
        , caption = "Macroalgae parameters and variables",bootstrap_options = "hover")%>%
   pack_rows(MAG_Groups[1],
            min(which(MAG_tab$Group == MAG_Groups[1])),
            max(which(MAG_tab$Group == MAG_Groups[1])),
            background = '#ebebeb'
            ,color="black") %>%
   pack_rows(MAG_Groups[2],
            min(which(MAG_tab$Group == MAG_Groups[2])),
            max(which(MAG_tab$Group == MAG_Groups[2])),
            background = '#ebebeb'
            ,color="black") %>%
   kable_styling(MAG_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>% 
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "10em" ,color="black",bold = F) %>%
  column_spec(2, width_min = "20em" ,color="black") %>%
  column_spec(3, width_min = "15em" ,color="black") %>%
  # column_spec(4, width_min = "15em" ,color="black") %>%
  row_spec(1:4, background = 'white')  %>%
 scroll_box(width = "45em", height = "25em",fixed_thead = FALSE)
```

</center>

<br>

The presence of macrophytes also affects the diffusivity at the sediment-water interface. Diffusion is calculated as normal in all layers, but the diffusivity in the top-most layer only is multiplied by the factor *swi_deff*. *swi_deff* is calculated according to equation \@ref(eq:MAG-10), where *TAUB* is the linked variable, and `taubsensitivity` and `smothbm` are CANDI-AED parameters, set for each zone. The effect is that when macrophytes are present, they smother the sediment and prevent mixing with the water column.\

```{=tex}
\begin{equation}
If  \  \   \ mag \gt smothbm 
\\
taub = 0.001
\\
(\#eq:MAG-8)
\\
\end{equation}
```
```{=tex}
\begin{equation}

If  \  \  \ mac \gt smothbm  
\\
taub = 0.001
\\
(\#eq:MAG-9)
\\
\end{equation}
```
```{=tex}
\begin{equation}

swi \text{_} deff = (1.0 + TAUB \times 100) \times taubsensitivity 
\\
(\#eq:MAG-10)
\
\\ 
\end{equation}
```
```{r PhysChem4, echo = FALSE, out.width='50%', class = "text-image",fig.show='hold',fig.align='center', fig.cap = "Schematic of macrophyte model: $O_2$ injection into the sediment by seagrass roots and smothering by filamentous algae."}
knitr::include_graphics("images/23-sediment_biogeochemistry/PhysChemInteractions/O2Inject-12.png")

```

### Implementation within the AED framework

While the first part of the model description was about the scientific processes, this second part describes how CANDI-AED functions as one module within the AED framework.

The AED modules are designed to fit together to make a simulation with a desired degree of complexity. Most water quality modules separate reaction processes, and the AED framework is used to combine them. For example, AED Oxygen and AED Organic Matter can be combined with a hydrodynamic model. CANDI-AED is one of the most complex of the modules, in part because it has both reaction and transport components within the module. The processes are not compartmentalised into separate models, and so many switches and parameters are used within CANDI-AED to set the degree of complexity and links to the other modules.

#### Sediment-water coupling

The sediment and hydrodynamic models are coupled at the sediment-water interface. The AED model setup has separate functions for coupling the variables of the bottom-most cell of a hydrodynamic model to the top-most layer of the sediment model:

-   flux of solid (particulate) material onto the sediment surface, mmol m^-2^ d^-1^
-   concentration of dissolved substances in the bottom water, mmol m^-3^
-   flux of dissolved substances from the top sediment layer to or from the water, mmol m^-2^ d^-1^

Where there is any mismatch between the sediment and water column variables, the factor of `part_sed_scale` can be used. (`part_sed_scale` can be found in the file `aed_sdg_vars.csv`, as outlined in more detail [below](#SedSedimentWaterInterfaceBoundary).) For example, influx of the water column variable $PON$ can be distributed between sediment variables $PON_R$ and $PON_L$. `part_sed_scale` is a real number between zero and one. For any sediment variables linked to the same water column flux, the `part_sed_scale` values should sum to 1.

<br>

```{r dev-pics2, echo = FALSE, out.width='45%', fig.align='center', fig.cap = "Schematic of sediment water coupling interactions"}

knitr::include_graphics("images/23-sediment_biogeochemistry/SWIBoundary/LinkedVars-08.png")

```

Depending on the nature of the host hydrodynamic model, several configurations can be implemented:

<center>

```{r 23-pic6, echo=FALSE, fig.cap="Spatial resolution options available through AED. a) Water column studies have traditionally assigned a flux to the sediment water interface without resolving the sediment chemical concentrations by depth, though they can be resolved laterally. b) The 0D water column is the method used in most sediment diagenesis studies, and use of multiple sediment zones is an option available within AED.", out.width = '95%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/image6.png")
```

</center>

Fluxes of dissolved species occur between the sediment and water column. They are calculated from the concentration gradient at the sediment-water interface according to Fick's Law:

<center>

```{=tex}
\begin{equation}
F = D_{0}\frac{\Delta C}{\Delta x} = \frac{D_{0}}{\delta} (C_{bw} - C_{1})
(\#eq:biogeochem1)
\end{equation}
```
</center>

where $D_{0}$ is the diffusivity, $\delta$ is the thickness of the diffusive boundary layer at the sediment water interface and defined as the length scale of the first sediment layer, $C_{bw}$ is the bottom water concentration and $C_{1}$ is the concentration in the top sediment layer.

At the bottom of the domain ($x$ = `xl`) the model can be specified to have a fixed-concentration (`deepbc_mode` = 0) such that the concentration at $xl = C_{Bot}$, or it can be specified to have a zero-derivative (`deepbc_mode` = 1) defined as $\frac{dC}{dx}=0$ at $x = xl$.

AED environmental variables can be accessed by CANDI-AED, and in most cases they do not feed back to the water column (see section on [environmental variables](#SedEnvironmentalVariables)). Parameter-variables can also be accessed by CANDI-AED, for example *taub* and *w~00~*. *taub* affects the diffusivity and the sediment-water interface. The deposition rate *w~00~* can be linked to the water column rate using `swi_link_variable`.

#### Resolving sediment zonation

The sediment model can be set up with multiple zones. The setup is the same as in the description above, and the zone boundares are not necessarily coincident with the grid structure of the water. Using zones is a practical compromise between computational efficiency and capturing spatial heterogeneity in sediment properties and their fluxes.

For a 3D hydrodynamic model, the sediment zones are adjacent areas of the bottom of the domain, for example, stretches of a river from upstream to downstream. A 3D diagram is given in Figure \@ref(fig:dev-Grid1).

```{r dev-Grid1, echo = FALSE, out.width='100%', class = "text-image", fig.cap = "Schematic depicting sediment zone numerical approach, for a 3D hydrodynamic mesh."}
knitr::include_graphics("images/23-sediment_biogeochemistry/GridStructure/3DZones-05.png")

```

The equivalent 1D hydrodyanmic model has the water column simulated in vertical layers. The sediment zones are assigned as slices of sediment according to depth. A 1D diagram is given in Figure \@ref(fig:dev-Grid2).

```{r dev-Grid2, echo = FALSE, out.width='100%', class = "text-image", fig.cap = "Schematic depicting sediment zone numerical approach, for a 1D layered hydrodyanmic model, such as a lake.."}
knitr::include_graphics("images/23-sediment_biogeochemistry/GridStructure/1DZones-05-04.png")

```

<br>

The fluxes and concentrations in the water cells above the sediment are averaged for linked variables (Figure \ref(fig:dev-Grid3)).

```{r dev-Grid3, echo = FALSE, out.width='60%', class = "text-image",fig.show='hold',fig.align='center', out.width='60%', fig.cap = "Schematic displaying how water cells are averaged when using sediment zones with a 3D hydrodynamic model."}

knitr::include_graphics("images/23-sediment_biogeochemistry/GridStructure/ZoneAveraging-04.png")

```

```{r dev-Grid4, echo = FALSE, out.width='60%', class = "text-image",fig.show='hold',fig.align='center', fig.cap = "Schematic displaying how water cells are averaged when using sediment zones with a 1D hydrodynamic model."}
knitr::include_graphics("images/23-sediment_biogeochemistry/GridStructure/ZoneAveraging-05.png")
```

#### Numerical solution {#SedNumericalSolution}

The time-dependent processes in CANDI-AED are solved by an ordinary differential equation solver DVODE. The details of the numerical solution are explained in detail by Boudreau (1997). DVODE was originally developed by Brown et al. (1989) and was used in the predecessors of this code by Boudreau (1996) and Luff et al. (2000).

CANDI-AED sets up the inputs for DVODE as follows:

-   The **boundary conditions** are updated at the sediment-water interface and bottom boundary (which may be connected to other modules)

<center>

```{r Numeric1, echo=FALSE, fig.cap="", out.width = '30%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/SWIBoundary/SWIboundary-03.png")
```

</center>

-   The **rates** are calculated, for example, as a product of concentration and rate constants, for example,\

    <center>

    ```{=tex}
    \begin{equation}
    R_{FeSOx} = k_{FeSOx} [O_2] [FeS]
    \end{equation}
    ```
    </center>

-   The **reactions** are calculated, as the sum of the losses or gains from the rates of consumption or production, for example,

<center>

```{=tex}
\begin{equation}
\frac {\delta O_2}{\delta t} = - R_{OM} - R_{NH4Ox} + R_{gpp} - R{rsp} + {\ \ } ... 
\end{equation}
```
</center>

-   The **transport** is calculated, as the sum of advection, diffusion and biological mixing

<center>

```{r Numeric2, echo=FALSE, fig.cap="", out.width = '25%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/AdvectionDispersion/AdvectionDispersion-14.png")
```

</center>

For each time step, one numerical matrix is produced, which contains an entry at each depth. In the model code, the matrix is given the simple name `y`. There is a corresponding matrix, named `ydot`, which contains all of the changes in the variables between time steps. The matrix `y` contains:

-   Dissolved variables
-   Solid variables
-   Other variables such as pH, pe, ubalchg and IAPs
-   Reaction sums
-   Diffusion
-   Advection

For example, for the reaction of $O_2$ and $FeS$, the matrix contains a depth entry for concentration of both species. The reaction between them is calculated according to the reaction $R_{FeS}$, which also has an entry at each depth. The $O_2$ is subject to transport processes for solutes, and the $FeS$ for solids. The concentration of the reactants are updated, along with the concentration of products such as $Fe^{2+}$ and $SO_4^{2-}$. This process is then followed for all other species.

<center>

```{r Numeric3, echo=FALSE, fig.cap="Representation of the matrix structure of the ordinary differential equation solver for and example reaction. With an entry at each depth, the ODE solver calculates the concentrations of reactants from the reactions and transport. ", out.width = '75%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/ProgramStructure/Matrix-05.png")
```

</center>

##### Solute transport equations: subroutine LiqMotion

A subroutine in the model code, named `LiqMotion`, calculates the solute transport equations. It ultimately produces `ydot` for solutes, which is used in the difference between time steps. `ydot` is the sum of diffusion, advection and irrigation, which are calculated separately. They are also calculated differently for each depth layer.

###### The top layer of the array

###### The top layer of the sediment

###### All intermediate layers

###### The bottom-most layer

##### Solid transport equations: subroutine SolMotion

###### The top layer of the array

###### The top layer of the sediment

###### All intermediate layers

###### The bottom-most layer

#### Module program structure

The sediment diagenesis model CANDI-AED can be used within the AED framework in various ways. This includes a) how the SDG module links to other modules in simulating water variables, and b) how the module operates within the simulated domain.

The general structure of the program is shown in Figure \@ref(fig:23-pic7). The program is firstly initialised (including spin-up days if desired), then loops through the kinetic and equilibrium reactions for each time step and writes the resulting concentrations and rates at each depth to an output file. The kinetic reactions are solved by the VODE program (Brown et al. 1989) and the equilibrium reactions by the Simplex program.

<center>

```{r 23-pic7, echo=FALSE, fig.cap="Candi program workflow. The program is firstly initialised, then loops through the kinetic and equilibrium reactions for each time step and writes the resulting concentrations and rates at each depth to an output file. The kinetic reactions are solved by the VODE program and the equilibrium reactions by the Simplex program.", out.width = '75%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/ProgramStructure/ProgramStructure-04.png")
```

</center>

### Optional module links

In the hierarchy of AED modules, CANDI-AED lies towards the bottom: CANDI-AED can be dependent on many other modules for boundary conditions, and few other modules are dependent on the results of CANDI-AED. CANDI-AED can be run effectively as a stand-alone model or it can be linked to the AED modules.

As outlined in the [section on the sediment-water interface](#SedSedimentWaterInterfaceBoundary), concentration or solid particle flux at the sediment-water interface can supply the top layer of sediment with the value from the bottom most water column cell. A simulation using CANDI-AED can use the following other modules:

-   **aed_sedflux**:
-   **aed_tracer**
-   **aed_noncohesive**
-   **aed_oxygen**: $O_2$ concentration can be linked via `OXY_oxy`
-   **aed_carbon**: $CH_4$ and $DIC$ concentration can be linked via `CAR_ch4` and `CAR_dic`
-   **aed_silica**
-   **aed_nitrogen**: $NH_4^+$ and $NO_3^-$ concentration can be linked via concentration via `NIT_amm` and `NIT_nit`
-   **aed_phosphorus**: $PO_4^{3-}$ concentration can be linked via `PHS_frp`
-   **aed_organic_matter**: several solid organic matter fluxes can be linked for carbon, nitrogen and phosphorus. including `OGM_poc_swi`, `OGM_pon_swi`, and `OGM_pop_swi`
-   **aed_phytoplankton**: the macroalgae settings in CANDI-AED can be linked to environmental variables calculated in the phytoplankton module
-   **aed_totals**:

The environmental variables can be used to set parameters and configurations in CANDI-AED, as outlined in the [section below](#SedEnvironmentalVariables).

### Feedbacks to the host model

As mentioned above, there are few other modules that are dependent on the results of CANDI-AED, unless the dissolved fluxes at the sediment-water interface are switched on:

-   **aed_sedflux**: the parameter `sedflux_model` controls the setup of zones and whether the flux is taken from a constant value or from CANDI-AED
-   **aed_tracer**
-   **aed_noncohesive**
-   **aed_oxygen**: $O_2$ flux can be linked via `SDF_Fsed_oxy`
-   **aed_carbon**: $CH_4$ and $DIC$ flux can be linked via `SDF_Fsed_ch4` and `SDF_Fsed_dic`
-   **aed_silica**
-   **aed_nitrogen**: $NH_4^+$ and $NO_3^-$ fluxes can be linked via `SDF_Fsed_amm` and `SDF_Fsed_nit`
-   **aed_phosphorus**: $PO_4^{3-}$ flux can be linked via `SDF_Fsed_frp`
-   **aed_organic_matter**: several dissolved organic matter fluxes can be linked for carbon, nitrogen and phosphorus, including `SDF_Fsed_doc`, `SDF_Fsed_don` and `SDF_Fsed_dop`
-   **aed_phytoplankton**:
-   **aed_totals**:

### Variable summary {#VariableSummary}

CANDI-AED contains the largest number of variables of any AED water quality module. This section outlines the state variables and diagnostic variables, including those that are compulsory and optional.

#### State variables

A list of the state variables is given in table \@ref(tab:VarsTable). The compulsory variables are key to the biogeochemical processes in CANDI-AED and there is a switch to stop the model running if they are not requested in the variables input csv. Three variables ($pH$, $pe$ and $ubalchg$) are not requested in the input file but they are always calculated and written to the output files. The non-compulsory variables are optional and if they are not included, the model will be allowed to run. However, the results of the calculations may be unpredictable and it is recommended that unless the user is familiar with this model, the variables be left in.

The organic matter variables are more straightforward in this regard, as the `OMModel` switch leaves the unused organic matter species unreacted. For example, if `OMModel` = 1, then $POML$ will be a reactive species and $POCL$ will not react. $POCL$ will still be transported and written to the output file, but its concentration will simply disperse from the initial condition.

<center>

```{r VarsTable, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
sediment_VarsTable <- read_xlsx("tables/23-sediment_biogeochemistry/BookdownSDGVars.xlsx", sheet="CompulsoryVars")
# View(sediment_VarsTable)
sediment_VarsTableGroups <- unique(sediment_VarsTable$Group)


  
  kable(sediment_VarsTable[,1:2],"html", escape = F, align = "c"
        , caption = "Full list of state variables.",bootstrap_options = "hover")%>%
  
    
  pack_rows(sediment_VarsTableGroups[1],
            min(which(sediment_VarsTable$Group == sediment_VarsTableGroups[1])),
            max(which(sediment_VarsTable$Group == sediment_VarsTableGroups[1])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(sediment_VarsTableGroups[2],
            min(which(sediment_VarsTable$Group == sediment_VarsTableGroups[2])),
            max(which(sediment_VarsTable$Group == sediment_VarsTableGroups[2])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(sediment_VarsTableGroups[3],
            min(which(sediment_VarsTable$Group == sediment_VarsTableGroups[3])),
            max(which(sediment_VarsTable$Group == sediment_VarsTableGroups[3])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(sediment_VarsTableGroups[4],
            min(which(sediment_VarsTable$Group == sediment_VarsTableGroups[4])),
            max(which(sediment_VarsTable$Group == sediment_VarsTableGroups[4])),
            background = '#ebebeb'
            ,color="black") %>%
   kable_styling(sediment_VarsTable, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>% 
    
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "10em" ,color="black",bold = F) %>%
  column_spec(2, width_min = "15em" ,color="black") %>%
  # column_spec(3, width_min = "15em" ,color="black") %>%
#  column_spec(4, width_min = "15em") %>%
  row_spec(1:4, background = 'white')  %>%
 scroll_box(width = "25em", height = "27em",fixed_thead = FALSE)
```

</center>

#### Diagnostic variables

There is no exhaustive list of the diagnostic variables in CANDI-AED. As described in the [outputs section](#SedOtherVariablesExtrasSed) below, most rates and factors can be written to an output .sed file using the list `morevariables`.

Another set of sheet variables is calculated within CANDI-AED and written to the netcdf output (table \@ref(tab:DiagVarsTable)). In order to set the amount of data written to the netcdf and prevent the creation of excessively large files, the parameter `diag_level` is used to set the number of diagnostics.

<center>

```{r DiagVarsTable, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
sediment_DiagVarsTable <- read_xlsx("tables/23-sediment_biogeochemistry/BookdownSDGVars.xlsx", sheet="Diag")
# View(sediment_VarsTable)
sediment_DiagVarsTableGroups <- unique(sediment_DiagVarsTable$Group)


  
  kable(sediment_DiagVarsTable[,1:4],"html", escape = F, align = "c"
        , caption = "List of diagnostic variables written to the netcdf.",bootstrap_options = "hover")%>%
  
    
  pack_rows(sediment_DiagVarsTableGroups[1],
            min(which(sediment_DiagVarsTable$Group == sediment_DiagVarsTableGroups[1])),
            max(which(sediment_DiagVarsTable$Group == sediment_DiagVarsTableGroups[1])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(sediment_DiagVarsTableGroups[2],
            min(which(sediment_DiagVarsTable$Group == sediment_DiagVarsTableGroups[2])),
            max(which(sediment_DiagVarsTable$Group == sediment_DiagVarsTableGroups[2])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(sediment_DiagVarsTableGroups[3],
            min(which(sediment_DiagVarsTable$Group == sediment_DiagVarsTableGroups[3])),
            max(which(sediment_DiagVarsTable$Group == sediment_DiagVarsTableGroups[3])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(sediment_DiagVarsTableGroups[4],
            min(which(sediment_DiagVarsTable$Group == sediment_DiagVarsTableGroups[4])),
            max(which(sediment_DiagVarsTable$Group == sediment_DiagVarsTableGroups[4])),
            background = '#ebebeb'
            ,color="black") %>%
   kable_styling(sediment_DiagVarsTable, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>% 
    
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "10em" ,color="black",bold = F) %>%
  column_spec(2, width_min = "15em" ,color="black") %>%
  # column_spec(3, width_min = "15em" ,color="black") %>%
#  column_spec(4, width_min = "15em") %>%
  row_spec(1:4, background = 'white')  %>%
 scroll_box(width = "40em", height = "27em",fixed_thead = FALSE)
```

</center>

#### Environmental variables {#SedEnvironmentalVariables}

There are several 'environmental' variables, that can be linked from the water column model or other AED modules, to set conditions in the sediment model (table \@ref(tab:EnvVarsTable)). For example, the set of $MPB$ and $MAG$ variables is used to set the $O_2$ penetration into the upper sediment. Another example is that water column salinity from a model such as TUFLOW-FV can be used to set the sediment-water interface concentration of salt at the sediment-water interface, which is used to calculate *F~Sal~*. Salinity is set as a permanently linked variable from the water column AED model, without using `water_link`.

Water column salinity is also used to set the boundary concentration of sulfate. Although most dissolved variables are simulated using mmol L^-1^, salinity uses practical salinity units (PSU). PSU has a standard concentration of solutes, primarily chloride, sodium, magnesium and sulfate. One PSU correlates to approximately 0.9 mmol L^-1^ sulfate, and so the CANDI-AED $SO_4^{2-}$ boundary is fixed at 0.9 $\times$ salinity.

<center>

```{r EnvVarsTable, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
sediment_EnvVarsTable <- read_xlsx("tables/23-sediment_biogeochemistry/BookdownSDGVars.xlsx", sheet="Enviro")
# sediment_EnvVarsTableGroups <- unique(sediment_VarsTable$Group)
  kable(sediment_EnvVarsTable[,1:4],"html", escape = F, align = "c"
        , caption = "List of environmental variables.",bootstrap_options = "hover")%>%
  
  kable_styling(MAG_tab, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>% 
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "10em" ,color="black",bold = F) %>%
  column_spec(2, width_min = "15em" ,color="black") %>%
  # column_spec(3, width_min = "15em" ,color="black") %>%
#  column_spec(4, width_min = "15em") %>%
  row_spec(1:4, background = 'white')  %>%
 scroll_box(width = "45em", height = "25em",fixed_thead = FALSE)
```

</center>

\

### Parameter and option summary {#SedParameterAndOptionSummary}

The preceding details were an explanation of the scientific aspects of this diagenesis model. Practical implementation is outlined below.

The model is set up via the name list text file `aed.nml` using the module keyword `&aed_seddiagenesis` and a model namelist block termed `&aed_sediment` configured to select the `sediment_model` as "Dynamic" or "Dynamic2D" option. Once this is selected, the model will search for the `&aed_sedcandi` parameter block.

In addition to adding the code block to `aed.nml`, users must also supply a valid AED sediment parameter database file (`aed_candi_pars.csv` is the default name) and a variable definition file (`aed_sdg_vars.csv` is the default name). If external boundaries are required, then a `swibc_file` and a `deepbc_file` also need to be supplied. The four key places for the model setup are summarised in table \@ref(tab:3dev-sedBouIni)), along with the file locations, and whether the settings are different for a simulation with more than one zone.

Users can create a standard file in the correct format from the online [**AED parameter database**](https://aed.see.uwa.edu.au/research/models/AED/aed_dbase/db_edit.php) by selecting from the available sediment types of interest, downloading via the **"Make CSV"** button, and then tailoring to the simulation being undertaken as required.

<center>

```{r 3dev-sedBouIni, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
  options(knitr.kable.NA = "")
sediment_BoundaryInitial <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet="BoundaryInitial")
  kable(sediment_BoundaryInitial,"html", escape = F, align = "c", caption = "File locations for major inputs of parameters, initial conditions and boundary conditions.",
      bootstrap_options = "hover") %>%
  kable_styling(sediment_BoundaryInitial, bootstrap_options = "hover",
                full_width = T, position = "left",
                font_size = 12) %>%
  row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
  column_spec(1, width_min = "10em" ,color="black",bold = T) %>%
  column_spec(2, width_min = "15em" ,color="black") %>%
  column_spec(3, width_min = "15em" ,color="black") %>%
#  column_spec(4, width_min = "15em") %>%
  row_spec(1:4, background = 'white') %>%
  scroll_box(width = "50em", height = "27em",fixed_thead = FALSE)
```

</center>

#### Module settings: &aed_sed_candi

The parameters set at the module level are outlined in table \@ref(tab:param-nml). These parameters are the same for all zones.

<center>

```{r param-nml, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)

param_nml <- read_xlsx("tables/23-sediment_biogeochemistry/BookdownParams.xlsx", sheet = "aed_sed_candi" )
  param_nmlGroups <- unique(param_nml$Group)
kable(param_nml[,1:3],"html", escape = F, align = "c", caption = "Parameters in &aed_sed_candi.",
     bootstrap_options = "hover") %>%

  
   kable_styling(param_nml, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
  pack_rows(param_nmlGroups[1],
            min(which(param_nml$Group == param_nmlGroups[1])),
            max(which(param_nml$Group == param_nmlGroups[1])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(param_nmlGroups[2],
            min(which(param_nml$Group == param_nmlGroups[2])),
            max(which(param_nml$Group == param_nmlGroups[2])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(param_nmlGroups[3],
            min(which(param_nml$Group == param_nmlGroups[3])),
            max(which(param_nml$Group == param_nmlGroups[3])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(param_nmlGroups[4],
            min(which(param_nml$Group == param_nmlGroups[4])),
            max(which(param_nml$Group == param_nmlGroups[4])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(param_nmlGroups[5],
            min(which(param_nml$Group == param_nmlGroups[5])),
            max(which(param_nml$Group == param_nmlGroups[5])),
            background = '#ebebeb'
            ,color="black") %>%
  pack_rows(param_nmlGroups[6],
            min(which(param_nml$Group == param_nmlGroups[6])),
            max(which(param_nml$Group == param_nmlGroups[6])),
            background = '#ebebeb'
            ,color="black") %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "15em",color="black") %>%
 column_spec(2, width_min = "15em",color="black") %>%
 column_spec(3, width_min = "10em",color="black") %>%
 row_spec(1:3, background = 'white') %>%
 scroll_box(width = "700px", height = "400px",
            fixed_thead = FALSE)
```

</center>

#### Zone-specific parameters

The majority of CANDI-AED parameters are specified for each zone, using the input parameter file. The default name is `aed_candi_pars.csv` but the user can change this in `aed.nml` `&aed_sedcandi`.

##### *Physical domain and numerical parameters* {.unnumbered}

The physical properties of the sediment are given in Table \@ref(tab:param-Mixing) and \@ref(tab:param-Layers). In general, the SDG module was set to run with a time-step of 2 hours, resolve a vertical sediment depth of 80 cm using 50 layers.

Some parameters are set as constant values for all zones, while others are set to be variable in each zone. The variability of bioturbation and irrigation along the length of the lagoon is set using the depth of the effect ($x_1$, $x_2$ and $x_{irrig}$) rather than the intensity ($DB_0$, $\alpha_0$). Upper porosity $\rho_0$ is held at 0.9 (i.e. 90% porewater) for all zones and the lower porosity $\rho_{00}$ is based on measured data at sites in each zone.

<center>

```{r param-Mixing, echo=FALSE, message=FALSE, warning=FALSE}
param_Mixing <- read_xlsx("tables/23-sediment_biogeochemistry/BookdownParams.xlsx", sheet = "Mixing" )
kable(param_Mixing,"html", escape = F, align = "c", caption = "Sediment physical properties, given in aed_candi_parameters.csv.",
     bootstrap_options = "hover") %>%
   kable_styling(param_Mixing, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "10em") %>%
 column_spec(2, width_min = "20em") %>%
 column_spec(3, width_min = "10em") %>%
 row_spec(1:5, background = 'white') %>%
 scroll_box(width = "600px", height = "380px",
            fixed_thead = FALSE)
```

\

```{r param-Layers, echo=FALSE, message=FALSE, warning=FALSE}
param_Layers <- read_xlsx("tables/23-sediment_biogeochemistry/BookdownParams.xlsx", sheet = "Layers" )
kable(param_Layers,"html", escape = F, align = "c", caption = "Sediment layer setup parameters and mixing parameters that might vary by zone, to set simulation with differences in sediment quality across space, given in aed_candi_parameters.csv.",
     bootstrap_options = "hover") %>%
   kable_styling(param_Layers, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "10em") %>%
 column_spec(2, width_min = "20em") %>%
 column_spec(3, width_min = "10em") %>%
 row_spec(1:5, background = 'white') %>%
 scroll_box(width = "600px", height = "350px",
            fixed_thead = FALSE)
```

</center>

\

##### *Organic matter parameters* {.unnumbered}

Many of the key processes in the sediment are driven by organic matter oxidation. Fresh organic matter supplied to the sediment surface fuels the bacterially-driven primary redox reactions that constantly shift the chemical equilibrium in the sediment. Organic matter metabolism also releases organic N and P to the water column. The parameters are given in Table \@ref(tab:param-OM).

The organic matter model option number 2 was used in this project. This option has both particulate and dissolved organic phases, and simulates organic C, N and P as separate state variables.

<center>

```{r param-OM, echo=FALSE, message=FALSE, warning=FALSE}
param_OM <- read_xlsx("tables/23-sediment_biogeochemistry/BookdownParams.xlsx", sheet = "OM" )
kable(param_OM ,"html", escape = F, align = "c", caption = "Organic matter oxidation parameters, given in aed_candi_parameters.csv.",
     bootstrap_options = "hover") %>%
   kable_styling(param_OM, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "10em") %>%
 column_spec(2, width_min = "20em") %>%
 column_spec(3, width_min = "10em") %>%
 row_spec(1:5, background = 'white') %>%
 scroll_box(width = "600px", height = "300px",
            fixed_thead = FALSE)
```

</center>

<br>

\

##### *Secondary redox reactions* {.unnumbered}

Secondary redox reactions are formed from the by-products of primary redox reactions. The nitrogen redox cycle parameters are given in Table \@ref(tab:param-Nredox).

<center>

```{r param-Nredox, echo=FALSE, message=FALSE, warning=FALSE}
param_Nredox <- read_xlsx("tables/23-sediment_biogeochemistry/BookdownParams.xlsx", sheet = "NRedox" )
kable(param_Nredox ,"html", escape = F, align = "c", caption = "Nitrogen redox cycle parameters, given in aed_candi_parameters.csv.",
     bootstrap_options = "hover") %>%
   kable_styling(param_Nredox, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "10em") %>%
 column_spec(2, width_min = "20em") %>%
 column_spec(3, width_min = "10em") %>%
 row_spec(1:5, background = 'white') %>%
 scroll_box(width = "600px", height = "275px",
            fixed_thead = FALSE)
```

</center>

\
\
\
\

Secondary redox parameters for other processes are given in Table \@ref(tab:param-SecondaryRedox).

<center>

```{r param-SecondaryRedox, echo=FALSE, message=FALSE, warning=FALSE}
param_SecondaryRedox <- read_xlsx("tables/23-sediment_biogeochemistry/BookdownParams.xlsx", sheet = "SecondaryRedox" )
kable(param_SecondaryRedox ,"html", escape = F, align = "c", caption = "Secondary redox parameters, given in aed_candi_parameters.csv.",      bootstrap_options = "hover") %>%
   kable_styling(param_SecondaryRedox, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "10em") %>%
 column_spec(2, width_min = "20em") %>%
 column_spec(3, width_min = "10em") %>%
 row_spec(1:5, background = 'white') %>%
 scroll_box(width = "600px", height = "325px",
            fixed_thead = FALSE)
```

</center>

<br>

\

##### *Geochemistry parameters* {.unnumbered}

Equilibrium and precipitation/dissolution constants are given in Table \@ref(tab:param-Geochem).

<center>

```{r param-Geochem, echo=FALSE, message=FALSE, warning=FALSE}
param_Geochem <- read_xlsx("tables/23-sediment_biogeochemistry/BookdownParams.xlsx", sheet = "Geochem" )
kable(param_Geochem ,"html", escape = F, align = "c", caption = "Geochemistry and precipitation/dissolution parameters, given in aed_candi_parameters.csv.",
      bootstrap_options = "hover") %>%
   kable_styling(param_Geochem, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "10em") %>%
 column_spec(2, width_min = "20em") %>%
 column_spec(3, width_min = "10em") %>%
 row_spec(1:5, background = 'white') %>%
 scroll_box(width = "600px", height = "310px",
            fixed_thead = FALSE)
```

</center>

<br>

#### Hidden parameters: aed_readcandi.F90 {.unnumbered}

A set of parameters that the user is unlikely to need to change has been removed from the model input files and set in the model code (table \@ref(tab:param-Hidden)).

<center>

```{r param-Hidden, echo=FALSE, message=FALSE, warning=FALSE}
param_Hidden <- read_xlsx("tables/23-sediment_biogeochemistry/BookdownParams.xlsx", sheet = "hidden" )
kable(param_Hidden ,"html", escape = F, align = "c", caption = "Other parameters, which are not currently listed in the parameter input file.",
      bootstrap_options = "hover") %>%
   kable_styling(param_Hidden, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "30em") %>%
 row_spec(1:5, background = 'white') %>%
 scroll_box(width = "500px", height = "310px",
            fixed_thead = FALSE)
```

## Setup & Configuration

The model is set up via the name list text file `aed.nml` using the module keyword `&aed_seddiagenesis` and a model namelist block termed `&aed_sediment` configured to select the `sediment_model` as "Dynamic" or "Dynamic2D" option. Once this is selected, the model will search for the `&aed_sedcandi` parameter block.

An example `aed.nml` configuration block for the `aed_seddiagenesis` module that includes multiple active sediment zones, plus microphytobenthos (MPB) and seagrass links, is shown below:

<br>

```{fortran, eval = FALSE}
!###############################################################################
! aed_seddiagenesis
!-------------------------------------------------------------------------------
&aed_sediment
   sediment_model    = 'DYNAMIC'    ! engages the CANDI-AED model
   mpb_link_variable = 'PHY_mpb'    ! use to set link for MPB
   mag_link_variable = 'MA2_mag'    ! use to set link for macroalgae inputs
   mac_link_variable = 'MAC_mac'    ! use to set link for macrophyte root inputs
   swi_link_variable = 'NCS_swi'    ! use to set link for variable w00
/
&aed_sed_candi
  !-- Time Settings --!
   spinup_days     =  90
   spinup_dt       =   0.25
   driver_dt       = 900     
   substep         =   8       
  !-- Zones details --!
   n_zones         = 5
   active_zones    = 104,103, 53, 23, 12,11
   zone_types      =   5, 3, 3, 1, 1
  !-- General Setup Options  --!
   dbase           = './AED/aed_candi_params.csv'
   vars_files      = './AED/aed_sdg_vars.csv'
   geochem_file    = './AED/aed_geochem_pars.dat'
  !-- Sediment Boundary Conditions  --!
   swibc_mode      = 0      ! previously ibc2
   deepbc_mode     = 1      ! previously ibbc
   swibc_file      = './AED/aed_sediment_swibc.dat'
   deepbc_file     = './AED/aed_sediment_deepbc.dat'
   swibc_filevars  = '' ! 'oxy', 'nit', 'amm', 'frp', 'poml'   ! from_bc_file
   deepbc_filevars = '' ! ,OXY_oxy,                            ! use_deep_bc
   flux_scale      = 1
  !-- Initial Conditions  --!
   SolidInitialUnit= 'mmolLsolid'
   OMInitMethodL   = 'LI_I'
   OM_topL         = 1
   OM_minL         = 0.9
   OM_cfL          = 0.6
   InitMinDepthL   = 99
   OMInitMethodR   = 'LI_I'
   OM_topR         = 1
   OM_minR         = 0.3
   OM_cfR          = 0.6
   InitMinDepthR   = 50
   POMVR           = 0.3
  !-- Outputs  --!
   diag_level         = 10
   output_profiles    =  .TRUE.
   morevariables      = 'Rgpp','Rrsp','FO2','f_an'
   output_diag_vars   = 'mpb','oxy','amm','docl','pocl','docr'
   n_ddpths           = 2
   output_diag_depths = 1.0, 5.0  ! cm below swi
/
```

<br>

<br>

```{block2, pars-text, type='rmdnote'}
In addition to adding the above code block to `aed.nml`, users must also supply a valid AED sediment parameter database file (`aed_candi_pars`) and a variable definition file ('aed_sdg_vars'). Both of these files file must be supplied in either `CSV` format. 
Users can create a standard file in the correct format from the online [**AED parameter database**](https://aed.see.uwa.edu.au/research/models/AED/aed_dbase/db_edit.php) by selecting from the available sediment types of interest, downloading via the **"Make CSV"** button, and then tailoring to the simulation being undertaken as required. Carefully check the parameter units and values!
```

<br>

<br>

The sections above described the scientific and technical processes within the model, the links to other AED modules, and the range of variables and parameters. This section outlines how the model treats the initial condition and boundary conditions.

### Setup details and examples

This section contains a set of steps for how to set up a simple simulation. Due to the large number of variables and parameters in CANDI-AED a full set of example parameters is not presented here. However, the lists of the tables above all contain example values that the user can consult.

#### Hydrodynamic model setup details and examples

The duration of the sediment model simulation is set within the parameters of the hydrodynamic driver model. For example, if using GLM, the duration is set in `glm.nml`, `&time`, then either if `timefmt` = 2, the `start` and `stop` times are set, or if `timefmt` = 3, the `num_days` parameter sets the number of days. If using TUFLOW-FV, the time is set in the .fvc file, where the `start time` and `stop time` are set as dates.

The sediment zone setup of the hydrodynamic model is also set in the hydrodynamic model parameters. The sediment model can be run without any zones, however, if zones are required, the user will need to make sure the zone numbers are aligned.

#### Module settings and examples

Recommendations for the module-level settings, and information about the necessary modules that are related to CANDI-AED are given below. Please see the [parameter section](#SedParameterAndOptionSummary) above for the full list of names, units and examples of the module-level parameters.

##### *&aed_sed_candi* {.unnumbered}

Most parameters in the &aed_sed_candi block can be left at the default values. The `substep` and `driver_dt` can be adjusted to optimise the simulation for speed and numerical stability. If the simulation is stable, then `spinup_days` are not required. The list of `active_zones` should correspond to the zones set in the hydrodynamic model and the parameter files. Users should check that the file paths and file names of the parameter and variable files are correctly set up in the user's directories. The boundary condition parameters and the initial condition parameters can be left at their default values. From the output settings, the `output_profiles` can be set to *TRUE* so that the .sed files are written. Any `morevariables` and `output_diag_vars` are set here.

##### *&aed_models* {.unnumbered}

The user should ensure that `aed_seddiagenesis` is included in the `models` list. It must come after `aed_totals`.

##### *&aed_sedflux* {.unnumbered}

The parameter `sedflux_model` must be set to `Dynamic2d` for the sediment-water column models to be coupled.

##### *&aed_sed_const2d* {.unnumbered}

If the sediment model is running in fewer zones than the water column model, then the water column model takes constant values for the sediment fluxes. The fluxes are listed in this code block.

##### *&aed_sediment* {.unnumbered}

Some of the names of the linking variables are listed here, for example, `mpb_link_variable` = `PHY_mpb`, `mag_link_variable` = `MA2_mag`, `mac_link_variable` = `MAC_mac` and `swi_link_variable` = `NCS_swi`.

##### *Fsed linking variables* {.unnumbered}

The names of Fsed linking variables are set in the various nml blocks. If a linking flux is desired, the user should check that each linking variable is set properly in each block. For example, to link the *diss_flux_link* for $O_2$, the parameter `Fsed_oxy_variable` is set to `SDF_Fsed_oxy` in `aed_oxygen` and also in the variables .csv file. This applies similarly for `&aed_nitogen` `Fsed_amm_variable` = `SDF_Fsed_amm` and `SDF_Fsed_nit`.

#### Zone-specific parameter details and examples

Recommendations for the zone-level settings are given below. All parameters can be set for each zone and therefore used to configure zones differently. If only one zone is needed, then the first set of parameters are used. Please see the [parameter section](#SedParameterAndOptionSummary) above for the full list of names, units and examples of the zone-level parameters.

##### *Physical domain and numerical parameters* {.unnumbered}

The parameters `xl` and `maxnpts` set the depth of the simulation and the number of layers. `xl` is often set to around 30 cm. This is deeper than biological mixing and it is around the depth of a field measurement such as a grab sample or core. Deeper than around 50 cm, most diagenetic processes are stable compared to the upper layers and so it is unlikely that a simulation deeper than 50 is necessary. The number of layers, `maxnpts` can be around 50. The simulation will require enough layers to resolve a profile of the upper sediment layers and few enough layers to keep the output file size to a minimum. Some tuning of `maxnpts` may be necessary if the simulation is unstable.

The porosity parameters `p0` and `p00` are key parameters that can make a simulation specific to a study site. For example, a sandy sediment can have low porosity whereas a clay sediment can have very high porosity. Similarly, the sedimentation rate `w00` has a large effect on the dynamics of a simulation and it is an important parameter for making a simulation specific to a study site. Unfortunately, it is rarely measured in the environment.

Bioturbation and bioirrigation depths (`x1`, `x2` and `xirrig`) have a large effect on the evolution of concentration profiles. If the water column is well oxygenated then these parameters will influence the oxic depth. These parameters also reflect a study site. If an environment has a healthy benthic environment, there will be a lot of mixing, whereas a eutrophic, hypersaline or anoxic environment will have little or no mixing.

##### *Organic matter parameters* {.unnumbered}

The `OMModel` parameter is used to set the complexity of the organic matter setup. The simplest is `OMModel` = 1, which can also be most easily compared with published studies. Once the `OMModel` is set, the oxidation rate parameters are set. For example, for `OMModel` 1, the oxidation rates of the labile and refractory phases can be set to approximately `poml2dic` = 100 and `pomr2dic` = 1 y^-1^. (A literature survey of this parameter can be found in Paraska et al. 2014.)

##### *Secondary redox parameters* {.unnumbered}

Default parameters can be used for nitrogen redox parameters, unless a simulation is specifically trying to tune to concentrations of nitrogen redox species. Other secondary redox parameters can be left at the default settings.

##### *Geochemistry parameters* {.unnumbered}

The key geochemistry paramter to select is `rxn_mode` (see section [*Mineral precipitation and ageing*](#SedMineralPrecipitationAndAgeing) and \@ref(tab:Geochem-table4) above). `rxn_mode` = 0 turns off precipitation, ageing and *pH* and is therefore the simplest setup. `rxn_mode` = 2 allows the constant precipitation and ageing of minerals, without any equilibrium.

#### Boundary setting details and examples

##### *Sediment-water interface boundary* {.unnumbered}

Boundary condition settings are described below in the [*Boundary settings* section](#SedBoundarySettings), including example fluxes and concentrations. The simplest boundary setting is *default_vals* rather than linked variables or the boundary file. In general, oxidised variables are present at the sediment-water interface, for example, $O_2$, $NO_3^-$, $MnO_2$ and $Fe(OH)_{3A}$. (Reduced species can be used in the initial profiles.) $SO_4^{2-}$ is generally at the sediment-water interface in a saline environment and not in a freshwater environment.

#### Initial condition setup details and examples

Initial condition settings are described below in the [*Sediment initialisation options* section](#SedInitialisation), including examples of concentrations (of solids or in porewater). The simplest setup can start with low or zero concentrations and allow the profiles to evolve. In general, reduced species will be present in the sediment, for example $NH_4^+$, $Mn^{2+}$, $Fe^{2+}$, $H_2S$ and $CH_4$. The exact concentrations can be set for each variable according to any available measured data.

### Sediment initialisation options {#SedInitialisation}

Before the first time step starts, CANDI-AED loads in an initial concentration profile for each state variable. There is also a period of 'negative time' before the full reactive transport system begins.

#### Initialisation of concentration profiles {.unnumbered}

Concentration profiles are set using one of five options (figure \@ref(fig:Init-5b)): 

 - constant, `I_CON`
 - exponential decrease, `I_EXP`
 - linear, concentrated at the sediment-water interface, `I_SWI`
 - linear, concentrated at depth, `I_DEP`
 - a spike, `I_SPK`

The user sets the initial concentration, a depth where it changes `InitMinDepth` and a proportion (0 to 1) of the initial concentration `_Min`. 

The constant option as a constant concentration from the sediment-water interface to the deepest layer. This is the default for most variables if no other parameters are set. 

The exponential option as the initial concentration at the sediment-water interface and an exponential decay to the bottom of the sediment. This is the default for most labile organic matter species.  

The first linear option has the initial concentration at the sediment-water interface and a linear decrease in concentration to the `InitMinDepth`, where the concentration is the `_Min` parameter multiplied by the initial concentration. This is the default for most refractory organic matter species, as wel as $FeOH_A$ andd $MnO_{2A}$. 

The other linear option has the initial concentration at the `InitMinDepth`, and the a linear increase in concentration to the sediment-water interface, where the concentration is the `_Min` parameter multiplied by the initial concentration. This is the default for $NH_4^+$ and $Fe^{2+}$. 

The spike option has a peak of concentration at the depth `InitMinDepth`, a linear increase from the sediment-water interface to that depth, and a linear decrease from that depth to 2 $\times$ the `InitMinDepth`, below which the concentration is zero. The concentration peak is the initial concentration, and the concentration at the sediment-water interface and the deep sediment is the `_Min` parameter multiplied by the initial concentration. This option is like a combination of `I_SWI` and `I_DEP`. Three examples of the spike initial profile are shown below, depending on the depth of `InitMinDepth` (\@ref(fig:Init-5c)). 

The concentration for each variable is set in the input file `aed_sdg_vars.csv`, in the column `initial_vals` (\@ref(fig:Init-6b)). If the simulation has several zones, then the initial concentration for the zone is set in other columns in the csv, under the column headings `initial_01`, `initial_02` etc. If there is no initial column in the csv for a particular zone, then the zone defaults to the value in the `initial_vals` column. The parameters for the depth and proportion are set in `aed.nml`. 



```{r Init-5b, echo=FALSE, fig.cap="Schematic of four of the options for initial profiles.",fig.show='hold',fig.align='center' ,out.width = '75%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Initialisation/InitialProfile-09.png")
```

\


```{r Init-5c, echo=FALSE, fig.cap="Schematic of three possible ways the spike initial profile can appear, depending on the depth parameter `InitMinDepth`.",fig.show='hold',fig.align='center' ,out.width = '100%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Initialisation/InitialProfile-08.png")
```

\




```{r Init-6b, echo=FALSE, fig.cap="Simplified image of the `aed_sdg_vars.csv` input file. For initialisation, the `variables` need to be set, along with default `initial_vals` or zones-specific `initial_01`, `initial_02` etc. columns.",fig.show='hold',fig.align='center' ,out.width = '100%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/aed_sdg_vars/aed_sdg_vars-04.png")
```

\

#### Initialisation of organic matter profiles

For organic matter, the initial profiles are set to either the linear or exponential options. For the linear option, the initial concentration is still taken from `initial_vals` in `aed_sdg_vars.csv`. The upper and lower concentrations are set as percentages of the initial concentration using the parameters `OM_top` and `OM_Min`, which are found in the `&aed_sed_candi` block in the `aed.nml` namelist. There are separate parameters for labile and refractory organic matter, so the parameter names are in fact `OM_topL` and `OM_topR`, and `OM_minL` and `OM_minR`. These parameters are real numbers as decimals between 0 and 1. The `OM_Top` fraction applies to the surface concentration. The `OM_Min` applies at the depth `InitMinDepth` (cm), which is also set in the `&aed_sed_candi` block.

For the exponential option, the initial concentration decreases according to equation \@ref(eq#init-7b):

```{=tex}
\begin{eqnarray}
concentration = initial vals \times e^ {- OM cfl \times depth} 
(\#eq:init-7b)
\end{eqnarray}
```
The shape of the curve is determined by the parameter `OM_cfl`, which is found in the `&aed_sed_candi` block.

#### Timing switches {#SedTimingSwitches}

The duration of the simulation is determined by the hydrodynamic model and its parameter inputs. Other timing parameters are set as outlined below.

##### *Substep* {.unnumbered}

The general task for the reactive transport model is to calculate a change in concentration over time, as a function of a change in concentration over space and as a result of chemical reactions. The general form of this is given in equation \@ref(eq:tim-1). When this is solved numerically, it is discretised by small increments of time, which are referred to as 'time steps' (Δt in equation \@ref(eq:tim-2)).

To allow the part of the model that solves these differential equations to process a large change in concentration over space or because of a large reaction, the timestep can be multiplied by a substep factor. Examples of where this is needed include a sharp concentration gradient over a small area, or in influx of a large amount of highly reactive substances.

```{=tex}
\begin{eqnarray}
 \frac {\delta C} {\delta t} = \frac {\delta C} {\delta x} + reactions
(\#eq:tim-1)
\\
\\
 {\Delta C} = (\frac {\Delta C} {\Delta x} + reactions) \times \Delta t
(\#eq:tim-2)
\end{eqnarray}
```
##### *Spinup substep* {.unnumbered}

In some simulations the initial timesteps have a high probability that the numerical model may crash while processing the initial condition (figure \@ref(fig:Timing-2)). A slow substep overcomes this problem, however, it has the disadvantage of producing a longer run time. The solution for this problem is the tool `firststeps` (y), which sets an initial time period during which a slower substep (`substep_0`) is used, and after which the normal substep (`substep_1`) is used. The switch `timeswitch` must be set to 1 for this feature to work, otherwise the normal `substep` applies.

```{r Timing-2, echo=FALSE, fig.cap="Schematic of the timestep across time. A different timestep (`substep_0`) can be used for a period (`firsteps`) until another timestep is used (`substep_1`).   ",fig.show='hold',fig.align='center', out.width = '75%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/TimeSwitches/TimeSwitches-02.png")
```

##### *Transport without reaction* {.unnumbered}

Another tool to help move through the spinup period is to allow transport processes to occur before the reactions start (figure \@ref(fig:Timing-3)). This evens out sharp chemical concentration gradients. This allows the simulation to start before the proper time of the simulation, represented in \@ref(fig:Timing-3) as negative time. The parameter `spinup_days` is the number of timestep integers of negative time, and `spinup_dt` is the length of each timestep (h). No fluxes at the sediment-water interface occur during negative time.

```{r Timing-3, echo=FALSE, fig.cap="Schematic of the timestep over time. The initial period up to `spinup_days` occurs in 'negative time', at a duration of `spinup_dt`.",fig.show='hold',fig.align='center', out.width = '75%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/TimeSwitches/TimeSwitches-03.png")
```

##### *Evolution of concentration profiles* {.unnumbered}

After initialisation, the concentration profile develops through the course of the simulation. Some of the most common ways for the profile to develop are given in figure \@ref(fig:Init-7b).

```{r Init-7b, echo=FALSE, fig.cap= "Examples of how the concentration profiles develop after initialisation. Left: the variable accumulates at the surface. Centre: a concentration spike develops. Right: the variable accumulates in the deep sediment.  ",fig.show='hold',fig.align='center' ,out.width = '100%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Initialisation/ProfileEvolution-05.png")
```

### Boundary settings {#SedBoundarySettings}

#### Sediment-water interface boundary {#SedSedimentWaterInterfaceBoundary .unnumbered}

##### *Static option:* default_vals {.unnumbered}

For the boundary at the sediment-water interface, the value is set using the `default_vals` column in the input file `aed_sdg_vars.csv` (figures \@ref(fig:SWI-1), \@ref(fig:SWI-2)). These values are constant in time for the entire simulation. Solids are a deposition flux (mmol m^-2^ y^-1^) and solutes are a concentration (mmol m^-3^). Variables may have other switches or options applied to them, as outlined below, but otherwise the values default to this. To enable this option, set `swibc_mode` = 1.

```{r SWI-1, echo=FALSE, fig.cap="Sediment-water interface boundary inputs that are constant in time are set with `default_vals`. ",fig.show='hold',fig.align='center' ,out.width = '60%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/SWIBoundary/default_vals-10.png")
```

\

```{r SWI-2, echo=FALSE, fig.cap="Simplified image of the `aed_sdg_vars.csv` input file. The `variables` and `default_vals` need to be set for each row of the csv.",fig.show='hold',fig.align='center' ,out.width = '100%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/aed_sdg_vars/aed_sdg_vars-02.png")
```

##### *Dynamic options:* linked variables {.unnumbered}

There are four methods for linking variables to the water column model: - particle deposition flux using `part_sed_link` - solute flux using `diss_flux_link` - bottom water concentration using `water_link` - fixed linked variables

These do not need to be set for any variable, but they may be set for any variable. `part_sed_link` and `water_link` must be set for solids and solutes and cannot both be assigned to one variable. If a variable is not set to be linked, the csv contains double inverted commas (" ") and the value is set to the `default_vals`. To enable this option, set `swibc_mode` = 0.

As mentioned above, where there is any mismatch between the sediment and water column variables, the factor of `part_sed_scale` can be used. `part_sed_scale` can be found in the file `aed_sdg_vars.csv`. `part_sed_scale` is a real number between zero and one. For any sediment variables linked to the same water column flux, the `part_sed_scale` values should sum to 1.

If a variable is linked, it is linked for all zones.

Environmental variables have a fixed link to the water column: $Salinity$, $Temperature$ and $\tau_b$. These variables are key components of hydrodynamic models but do not undergo any rates of change in the sediment model. They are used to affect other sediment processes. As mentioned in the [*Environmental variables* section](#SedEnvironmentalVariables), the $SO_4^{2-}$ concentration (mmol L^-1^) is fixed at 0.9 $\times$ salinity (PSU).

```{r SWI-3, echo=FALSE, fig.cap="Schematic of possible linked variable options at the sediment-water interface. Particle fluxes and bottom water concentraitons go from the water column to the sediment. Sediment solutes can also be set to flux to the water column.",fig.show='hold',fig.align='center' ,out.width = '60%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/SWIBoundary/LinkedVars-08.png")
```

```{r SWI-4, echo=FALSE, fig.cap="Simplified image of the `aed_sdg_vars.csv` input file. Linking variables may be set in the same row as the variable, or left blank with inverted commas.",fig.show='hold',fig.align='center' ,out.width = '100%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/aed_sdg_vars/aed_sdg_vars-05.png")
```

##### *Dynamic options:* swibc file {#DynamicOptionsSwibcFile .unnumbered}

A dynamic boundary can be set for any variable and any zone using the text file swibc.dat (\@ref(fig:SWI-5)). For a variable to be activated to take its boundary from swibc.dat, it must be included in the list `swibc_filevars`, which is a parameter in the `&aed_sed_candi` block. The specific filename and filepath are set in the parameter `swibc_file` in the `&aed_sed_candi` block. The swibc file has the following features:

-   **timestep** column: The number of the point where anything changes from the previous entry.\\
-   **time** column: The time (days) at which a change occurs. Unlike a time series, each time step does not need to be listed. However, if there is a change in any one of the variables, an entry must be included in that row for all of the variables.\\
-   **forzone** row: The zone for which the variable in this file applies. If a zone is not listed, then the variable in that zone will take the `default_vals`.\\
-   **variable** columns: The value at the boundary until the next point of change.\\
-   **w00h00** columns: if any column is labelled `w00h00`, then a dynamic poreflux is applied to that zone.\\

The benefit of using this file is that a time-changing boundary can be applied for a selection of variables and zones, while leaving other combinations of variables and zones unchanged. To enable this option, set `swibc_mode` = 10.

```{r SWI-5, echo=FALSE, fig.cap="Schematic of the sediment-water interface boundary supplied by the `swibc.dat` text file.",fig.show='hold',fig.align='center' ,out.width = '60%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/SWIBoundary/swibcdat-09.png")
```

```{r SWI-6, echo=FALSE, fig.cap="Simplified image of the swibc.dat text file.",fig.show='hold',fig.align='center' ,out.width = '85%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/SWIBoundary/swibcdat-01.png")
```

##### *Special organic matter flux* {.unnumbered}

The sediment model has a unique boundary condition, where a large influx of organic matter occurs for a limited time (figure \@ref(fig:Timing-1)). This represents, for example, the operation of a fish farm, a large algal bloom or an illegal dump of organic waste into a lake. This is a convenient tool for looking at the effect of large amounts of organic matter on other chemical species, and it operates separately to the constant flux boundaries, the swibc input file, or the variables linked to the pelagic model. The flux defaults to 0.0 before the time parameter `fluxon` and after the time parameter `fluxoff` (y). The state variable switched on by this flux is $POM_{Special}$. It is not named as, for example, labile or refractory because the user can tailor its properties to the simulation. The breakdown parameter for $POM_{Special}$ is `pomspecial2dic` (y^-1^). The parameter ibc2 must be set to 10 or 0 for this feature to operate. The magnitude of $POM_{Special}$ flux is set in the file aed_sdg_vars.csv, in the default_vals column, and the default flux before `fluxon` and after `fluxoff` times is zero. To enable this option, set `swibc_mode` = 2.

```{r Timing-1, echo=FALSE, fig.cap="Schematic of $POM_{Special}$ flux over time. $POM_{Special}$ flux begins at time `fluxon` and stops at time `fluxoff`. The rate of the flux is `pomspecial`.",fig.show='hold',fig.align='center', out.width = '75%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/TimeSwitches/TimeSwitches-01.png")
```

##### *Combined organic matter flux and timestep switch* {.unnumbered}

A large influx of organic matter can have the same effect for the model as processing the initial condition, and so the large flux can crash the model well after the initial spinup (figure \@ref(fig:Timing-4)). Therefore, further substep parameters are available for slowing the model while it processes the organic matter. The parameter `substep_2` is used as the substep after `fluxon`. The parameter `substep_3` is used as an intermediate substep after `fluxoff`, until the time `justwaitabit`, to allow the final effects of the organic matter influx to be processed, but allowing the possibility of running faster if not at the speed of `substep_1`. Further *substep* options are explained in the [timing switches section](#SedTimingSwitches).

```{r Timing-4, echo=FALSE, fig.cap="Schematic of $POM_{Special}$ flux and timestep changes over time. Options are available to optimise a short run total time with a slow substep centred around organic matter influx.",fig.show='hold',fig.align='center', out.width = '75%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/TimeSwitches/TimeSwitches-04.png")
```

#### Bottom boundary {.unnumbered}

##### *Static option:* deep_vals {.unnumbered}

A boundary can be supplied beneath the deepest cell using `deep_vals` (\@ref(fig:Deep-1)). The value is taken from `deep_vals`, which is in the input file `aed_sdg_vars.csv` (\@ref(fig:Deep-2)). (`deep_vals` is the lower boundary equivalent of the sediment-water interface `default_vals`.) In most instances a bottom boundary would not be supplied, however, if there were a simulation with, for example, groundwater infiltration, then this option is used. If the value in `deep_vals` is negative, no boundary is applied, and so the file has a default -999 for all variables. To enable this option, set `deepbc_mode` = 1.

```{r Deep-1, echo=FALSE, fig.cap="Bottom boundary inputs that are constant in time are set with `deep_vals`. ",fig.show='hold',fig.align='center', out.width = '60%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/DeepBoundary/deep_vals-12.png")
```

```{r Deep-2, echo=FALSE, fig.cap="Simplified image of the `aed_sdg_vars.csv` input file. deep_vals are set if a value is present, greater than zero, corresponding to a variable.",fig.show='hold',fig.align='center' ,out.width = '100%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/aed_sdg_vars/aed_sdg_vars-03.png")
```

\

##### *Dynamic option:* deep_bc file {#DynamicOptionsDeepbcFile .unnumbered}

A dynamic boundary can be set for any variable and any zone using the text file deepbc.dat (\@ref(fig:Deep-3)). For a variable to be activated to take its boundary from deepbc.dat, it must be included in the list `deebc_filevars`, which is a parameter in the `&aed_sed_candi` block. The specific filename and filepath are set in the parameter `deepbc_file` in the `&aed_sed_candi` block. The deepbc file is the bottom boundary equivalent of the swibc file. It has the following features:

-   **timestep** column: The number of the point where anything changes from the previous entry.\\
-   **time** column: The time (days) at which a change occurs. Unlike a time series, each time step does not need to be listed. However, if there is a change in any one of the variables, an entry must be included in that row for all of the variables.\\
-   **forzone** row: The zone for which the variable in this file applies. If a zone is not listed, then the variable in that zone will take the `default_vals`.\\
-   **variable** columns: The value at the boundary until the next point of change.\\
-   **pf** columns: if any column is labelled `pf`, then a dynamic poreflux is applied to that zone.

The benefit of using this file is that a time-changing boundary can be applied for a selection of variables and zones, while leaving other combinations of variables and zones unchanged. To enable this option, set `deepbc_mode` = 10.

```{r Deep-3, echo=FALSE, fig.cap="Schematic of the bottom boundary supplied by the `deepbc.dat` text file.",fig.show='hold',fig.align='center', out.width = '60%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/DeepBoundary/DeepBC-11.png")
```

## Output & Post-processing

CANDI-AED has a range of output files and output file options.\

### Concentration of state variables: '.sed' files

The model writes output files for each state variable (chemical species, $IAP$s, $pH$, $pe$ and $ubalchg$) as text files, with a '.sed' extension. The files are written to the results directory, under the `candi_aed` subdirectory, and then in a separate subdirectory for each zone (with five significant figures), for example `.../results/candi_aed/00001/oxy.sed`. The '.sed' file has three rows of headers, then a row that lists the depths (cm). The first column is the sequence of output times (years). The other text entries are a matrix of the concentration values at each depth and each time (mmol L^-1^) (figure \@ref(fig:Output-1)). These can be read in using a data processing program such as R, Matlab or Excel for analysis and plotting \@ref(fig:Output-4).

```{r Output-1, echo=FALSE, fig.cap="Simplified version of the .sed output file.",fig.show='hold',fig.align='center', out.width = '85%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Outputs/Outputs-01.png")
```

```{r Output-4, echo=FALSE, fig.cap="Examples of the figures that can be plotted with the .sed file. Left: concentration-depth at one time step; centre: concentration-time at one depth; right: concentration-depth-time ",fig.show='hold',fig.align='center', out.width = '75%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Outputs/Outputs-04.png")
```

### Flux of state variables: 'swi_fluxes.sed'

The model writes output files for sediment-water interface fluxes, named `swi_fluxes.sed`. As with the concentration '.sed' files, the 'swi_fluxes' files are written to the results directory, under the `candi_aed` subdirectory, and then in a separate subdirectory for each zone (with five significant figures), for example `.../results/candi_aed/00001/swi_fluxes.sed`. The file has one row of the names of all state variables, and the second row for whether the variables are solids or solutes. The first column is the sequence of output times (years). Only the post-spinup times (positive time) are written to the file, as there are no fluxes during spinup. The other text entries are the flux values (mmol ^-2^ y^-1^), where a positive value is a flux into the sediment, and a negative value is a flux to the bottom water (figure \@ref(fig:Output-2)). As with the '.sed' files, these files can be processed and plotted (\@ref(fig:Output-5)).

```{r Output-2, echo=FALSE, fig.cap="Simplified version of the `swi_fluxes.sed` file.",fig.show='hold',fig.align='center', out.width = '85%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Outputs/Outputs-02.png")
```

```{r Output-5, echo=FALSE, fig.cap="Example of a figure that can be plotted from fluxes.sed ",fig.show='hold',fig.align='center', out.width = '75%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Outputs/Outputs-05.png")
```

### Other variables: '/Extras/ .sed' {#SedOtherVariablesExtrasSed}

Any variable listed in the `&aed_sed_candi` block parameter `morevariables` can be written if to a file if desired. Extra '.sed' files can be written to a subdirectory of the zone directory, for example `.../results/candi_aed/00001/Extras/RO2.sed`. The file has the same format as the other '.sed' files, however, it does not record values during the spinup period (figure \@ref(fig:Output-3)). The unit of the variable is dependent on what is used in the model code and there is no unit conversion before the file is written.

```{r Output-3, echo=FALSE, fig.cap="Simplified version of the extra .sed file. ",fig.show='hold',fig.align='center', out.width = '85%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/Outputs/Outputs-03.png")
```

### Depth profile settings: 'Depths.sed'

The model writes a text file that reports the initial depth profile settings, with one entry per layer (`maxnpts`). One file is written for each sediment zone and it is not updated throughout the simulation, for example `.../results/candi_aed/00001/Depths.sed`. The file can be useful in post-processing, including performing unit conversions. The outputs include

-   *Depths*: The depths of the grid layers from the sediment-water interface (cm)\
-   *Porosity*: The volume of water per volume of space, between 0 and 1\
-   *Irrigation*: The initial irrigation\
-   *Bioturbation*: The initial bioturbation\
-   *Layer_size_cm*: The sizes of the layers (cm)\
-   *Porewater_m3_per_m2*: The volume of porewater in each layer (m^3^ m^-2^)\
-   *Solids_m3_per_m2*: The volume of sediment solids (m^3^ m^-2^)\

### Mass balance calculation

The function *Mass* within the CANDI-AED code calculates the mass balance of some species.

Firstly it calculates *FOM_in*, which is the influx of

<!-- ### netcdf outputs: .nc{-} -->

<!-- As with other modules, the NetCDF files include the model outputs.  -->

<!-- ### Diagnostics: {-} -->

<!--   INTEGER  :: diag_level = 10                ! 0 = no diagnostic outputs\ -->

<!--                                              ! 1 = basic diagnostic outputs\ -->

<!--                                              ! 2 = flux rates, and supporitng\ -->

<!--                                              ! 3 = other metrics\ -->

<!--                                              !10 = all debug & checking outputs\ -->

<!-- ### Rates: rates.sed{-} -->

<!-- The model outputs the file `rates.sed` for each variable at every depth for the final time step. -->

## Case Studies & Examples

Several case studies are outlined below.

### Numerical mass balance and flux confirmation {-}

#### Mass check approach {-}

It is important to check that a simulation is conserving mass in the upper layers of the sediment because with the wrong settings, numerical errors can lead to the wrong calculation of sediment mass. One method for checking for errors involved a non-reactive solid and a non-reactive solute being introduced at the sediment-water interface (SWI) as 'tracers'. As they enter at the SWI, the tracers' fluxes and concentrations were converted to mass units, per 1 m$^2$. The concentrations in the sediment were similarly summmed per 1 m$^2$, across all depths. 

#### Mass check methods {-}

#### General settings {-}

The sediment model was run with no variables except one solid (organic matter) and one solute (N2) ‘tracer’, the reaction rates of which were set to zero. The tracers were applied as fluxes and concentrations at the sediment-water interface boundary. The simulation had 

-	no influx for 0.5 years to check that the model had equilibrated, then 
-	the boundary flux or concentration was applied for 1 year, then 
-	the boundary was switched off and the simulation run for a further year. 

The mass of influxing matter was calculated. Key parameters were changed one by one and a range of simulations was run, in order to check for possible mass balance errors. 

##### Settings investigated for errors {-}

Key parameters that were varied included 

-	either a fixed grid or exponentially increasing layers
-	grid spacing and time steps
-	solid sedimentation rates (ω00)
-	porosity (ρ0 and ρ00)
-	bioturbation depths and functions

##### SWI specific settings{-}

'aed.nml' settings:

- 'swibc_mode' = 10 (integer switch). The boundary to the file listed in 'swibc_file'.
- 'swibc_file' = 'swibc.dat' (filename). The sediment-water interface boundary condition (SWIBC) file in the aed_sdg folder.
- 'swibc_filevars' = 'pomspecial', 'n2' (variable names). The variables taken from the swibc file. 

'aed_candi_params.csv' settings:

- 'OMModel' = 1 (integer switch). The simplest organic matter reaction model. 
- 'pomspecial2dic' = 0 (y$^{-1}$). The 'pomspecial' that fluxes in is unreactive.

The 'pomspecial' flux was 1.0E-08 (mmol m$^{-2}$ y$^{-1}$) and the 'n2' concentration was 1.0E-08 (mmol m$^{-3}$) for most of the simulation. The model converted the concentration from mmol m$^{-3}$ to mmol L$^{-1}$.

At time 183 days (or 0.5 years) the flux and concentration increased to 1.0E+02 (mmol m$^{-2}$ y$^{-1}$) and 1.0E+05 (mmol m$^{-3}$). At 548 days (or 1.5 years) the flux and concentration returned to the very low number. 

##### Grid settings {-}

'aed_candi_params.csv' settings:

- 'job' = 0 (integer switch). An evenly-spaced grid.
- 'p0'  = 0.5 (fraction between 0 and 1). The porosity at the sediment-water interface.
- 'p00' = 0.5 (fraction between 0 and 1). The porosity at depth. These porosities were set to be the same, which gave the most reliable results.
- 'imix' = 2 (integer switch). The bioturbation set to have exponential decay rather than the sharper cutoffs with 'imix' = 1.
- 'maxnpts' = 41 (layers). With 'job' = 0 this must a multiple of the depth (cm) + 1 layer (for the zeroth layer), for example, 81, 201 or 401.
- 'xl' = 20 (cm). This depth was adequate for this simulation.
- 'xirrig' = 3 (cm). The depth of irrigation.
- 'xs' = 5. The shape of the decay curve for bioturbation. 

##### Mass calculation method {-} 

This section outlines the method for calculating the mass. There are many steps in the path to convert all of the quantities, starting with input parameters, calculations in the model, outputs from the model, and post-processing in R. An error in the assumptions behind or calculation of any of these steps would give a false idea of whether mass was being maintained. 

###### Parameter inputs {-}

Porosity is assigned in aed_candi_params.csv, using the parameters 'p0' and 'p00', setting porosity at the highest and lowest layers. 

###### Internal calculations {-}

'aed_secandi' calculates the variable 'poros' as an exponential decay from 'p0' to 'p00'. 

'poros' is used to calculate 'ps' (1 - poros). 'ps' and 'poros' are used to calculate the 'psp' and 'pps' as 'ps'/'poros' and 'poros'/'ps'. These are ratios of porosity that are used to convert solid and solute units elsewhere in candi. For example, as 'Fe(OH)3A' (solid) reduces to 'Fe$^{2+}$' (solute) it is multiplied by 'psp' to convert the units from mmol L $^{-1}$ solids to  mmol L $^{-1}$ (porewater). 

Two other variables are calculated within 'aed_sedcandi', using 'dh', which is the distance between sediment layers (cm). The porewater volume variable, 'porevol', is calculated as 'dh' / 100 $\times$ 'poros'. The solid volume variable, 'solvol', is calculated as 'dh' / 100 $\times$ 1 - 'poros'. These are then expressed as m$^3$ of porewater or solids per 1 m$^2$. 

###### Model outputs {-}

The model produces solid variables in units of mmol / L $^{-1}$ (solids) and solute variables in units of mmol L $^{-1}$ (porewater). These are written as concentrations over time, across all depths. 

The model also produces a time series of SWI fluxes. These are in units of mmol m$^{-2}$ y^{-1}, where the area component is total area of space, rather than the area of solids or porewater.

It also outputs 'porevol' and 'solvol', with entries per depth, but without time, since they do not change over time. 

###### R post-processing {-}

In processing the SWI fluxes file, the fluxes are converted via the following steps:

- raw output is in mmol m$^-2$ y$^-1$ / 1000
- mol m$^-2$ y$^-1$ $\times$ molecular mass 
- cumulative [t] g m$^-2$ = cumulative[t-1] g m$^-2$ + g m$^-2$ y$^-1$ [t] $\times$ (time[t] - time[t-1])

This gives a time series of influxed mass per 1 m$^-2$. 

The mass in the sediment is converted via the following steps:

- raw output is in mmol L$^-1$
- if the variable is a solid, the concentration is multiplied by solvol and 1000. 
  mmol L^-1 $\times$ m$^3$ of solids per 1 m $^2$
  = mmol solids per 1 m$^2$

  - if the variable is a solute, the concentration is multiplied by porevol and 1000. 
  mmol L^-1 $\times$ m$^3$ of porewater per 1 m$^2$
  = mmol porewater per 1 m$^2$

- mmol m$^{-2}$ / 1000 = mol m$^{-2}$ 
- mol m$^{-2}$ $\times$ molecular mass = g  m$^{-2}$
- this matrix of mass m$^{-2}$ is then summed over all depths at each time. This gives a similar time series of mass per 1 m$^{-2}$.
- for plotting, this is then carved into three sections, 0 cm to 1 cm, 1 cm to 2 cm, and 2 cm to the bottom. This does not affect the total mass. 

#### Mass check results {-}

##### Mass check base case  {-}

Simulations were run with a non-reactive solute and solid to check that the total influxing mass was equal to the mass across all layers of the sediment. Parameters and settings were changed until the masses were the same and the model could be used with confidence. The bottom left panel in each figure shows the sediment mass as the shaded area, and the influxed mass as the red line. An equal grid spacing was used to minimise possible large concentration differences between layers. The model was run quickly with 5-daily time steps, which had no effect on model error. The most important setting was that the initial bioturbation profile was set to the exponential decay option rather than the two-layer option (imix = 2 rather than 1, see [chapter 14.3.2.6](https://aquaticecodynamics.github.io/aed-science/sediment-biogeochemistry.html#model-description-8) for details) because it was found that the sharp cutoff of the two-layer option was prone to small numerical errors, when the concentration differences between layers were large. With the exponential decay setting, bioturbation decreased from its maximum value at the sediment water interface to 10 cm deep. 

```{r dev-masscheck1, echo = FALSE, out.width='75%', class = "text-image",fig.show='hold',fig.align='center', fig.cap = "Results of a simulation for $N_2$ as a tracer. Reaction rates were turned to zero and so concentration and mass changes are the result of transport only. After 6 months of spin-up, both species were introduced at the sediment-water interface for one year. The bottom left figure showed that the total sediment mass and influxing mass were the same."}

knitr::include_graphics("images/23-sediment_biogeochemistry/MassCheck/6P_n2_.png")

```


```{r dev-masscheck2, echo = FALSE, out.width='75%', class = "text-image",fig.show='hold',fig.align='center', fig.cap = "Results of a simulation for organic matter. Reaction rates were turned to zero and so concentration and mass changes are the result of transport only. After 6 months of spin-up, both species were introduced at the sediment-water interface for one year. The bottom left figure showed that the total sediment mass and influxing mass were the same."}

knitr::include_graphics("images/23-sediment_biogeochemistry/MassCheck/6P_pomspecial_.png")

```

##### Deliberate errors - high porosity {-}

Parameters 'p0' and 'p00' were set to 0.75. This resulted in a good solid balance but a higher sediment mass than influxed mass for solutes.

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/Bioturbation_hp.png" width=33%>
</p>

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/Porosity_hp.png" width=33%>
</p>

With high porosity, the solid mass reconciled well but the solute mass was higher in the sediment than from the influx. The maximum mass of solute was around 350 g m$^-2$. 

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/1P_Mass_n2_hp.png" width=67%>
</p>

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/1P_Mass_pomspecial_hp.png" width=67%>
</p>

##### Deliberate errors - low porosity {-}

Parameters 'p0' and 'p00' were set to 0.75. This resulted in a good solid balance but a higher sediment mass than influxed mass for solutes.

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/Bioturbation_lp.png" width=33%>
</p>

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/Porosity_lp.png" width=33%>
</p>

With low porosity, as with high porosity, the solid mass reconciled well but the solute mass was higher in the sediment than from the influx. The absolute mass of solid was the same with high or low porosity, however, the maximum mass of solute was approximately 100 g m$^2$. 

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/1P_Mass_n2_lp.png" width=67%>
</p>

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/1P_Mass_pomspecial_lp.png" width=67%>
</p>


##### Deliberate errors - porosity gradient {-}

Parameters 'p0' was set to 0.9 and 'p00' to 0.5. 

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/Bioturbation_pg.png" width=33%>
</p>

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/Porosity_pg.png" width=33%>
</p>

This created the clearest error. Reported sediment mass was much lower than influxed mass. 

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/1P_Mass_n2_pg.png" width=67%>
</p>

<p align="center">
<img src = "images/23-sediment_biogeochemistry/MassCheck/1P_Mass_pomspecial_pg.png" width=67%>
</p>

##### Bioturbation settings {-}

In this setup, bioturbation does not cause a problem with mass. However, future simulations with different bioturbation settings should be checked for mass conservation. One of the primary risks is having a sharp cutoff for bioturation at two depths.

The bioturbation intensity, 'DB0', might also affect the mass conservation. 


### Case study: Model benchmark example

#### Model setup

##### *Base simulation* {.unnumbered}

The model CANDI-AED was tested against the published modelling results of Van Cappellen and Wang (1996) in order to confirm that the model worked as expected, and also to demonstrate the flexibility of the model system. For this simulation, all boundary fluxes, bottom water concentrations, and rate constants were set to the same values as those in Van Cappellen and Wang (1996). The organic matter redox approach was set to Approach 2, by setting `OMApproach` = 2. The organic matter oxidation model was parameterized with no organic matter influx and with an oxidation rate that changes with depth (`OMModel` = 1, `VCW` = *TRUE*). The surface oxidation rate and depth attenuation were the same as Van Cappellen and Wang (1996). The irrigation coefficient was also set to be the same as in Van Cappellen and Wang (1996), using `irrigswitch` = 2. The mineral precipitation reactions were implemented only for MnCO3, FeCO3 and FeS, as per the equations in Van Cappellen and Wang (1996), rather than the larger set of precipitation reactions (`rxn_mode` = 4). Additionally, the ageing reactions of iron and manganese minerals were disabled. Ammonium adsorption was included and was the same as in Van Cappellen and Wang (1996), however, iron and manganese adsorption has not yet been included. An equivalent evenly spaced grid was used (`job` = 1), with 400 layers (`mxnpts`) to a depth of 20 cm (`xl`).

##### *Parameter variation* {.unnumbered}

The model was run with random variation in its parameter values, to assess the sensitivity of the outputs to the parameter values. The base simulation above was run again with variations of the uncertain parameters listed in Table 6 of Van Cappellen and Wang (1996). Rather than using the ranges that Van Cappellen and Wang used, the range was set using their chosen parameter value as the median. 100 simulations were run using a Matlab script, with individual parameter values automatically generated using a Monte Carlo-based method. The results of each simulation were then collated and the total results were assessed for distribution about the median value.

<center>

```{r param-VCW1, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
param_VCW1 <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet = "VCW1" )
kable(param_VCW1 ,"html", escape = F, align = "c", caption = "Parameters varied in the uncertainty analysis, within the ranges chosen by Van Cappellen and Wang (1996).",
      bootstrap_options = "hover") %>%
   kable_styling(param_Hidden, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "5em") %>%
 column_spec(2, width_min = "15em") %>%
 column_spec(3, width_min = "10em") %>%
 column_spec(4, width_min = "5em") %>%
 row_spec(1:5, background = 'white', color = "black") %>%
 scroll_box(width = "45em", height = "25em",fixed_thead = FALSE)
```

</center>

#### Results

##### *Concentrations* {.unnumbered}

The $O_2$ concentration profile matched the Van Cappellen and Wang simulation very closely, apart from the concentration at the deepest point, which was slightly higher than the field data (Figure \@ref(fig:VCW-1) (a)). The results for $O_2$ were very consistent about the median value. The slightly higher concentration of oxygen at the deepest point in this simulation may have carried through to inhibit denitrification, causing the concentration of $NO_3^-$ to be higher than in the Van Cappellen and Wang simulation (Figure \@ref(fig:VCW-1) (a)). $NO_3^-$ was similarly consistent about the median. The $NH_4^+$ field data points were lower than the median of these simulations but within the range of these simulations (Figure \@ref(fig:VCW-1) (c)). $pH$ was between 7 and 8, and both $NH_4^+$ and $pH$ had increasing ranges with depth (Figure \@ref(fig:VCW-1) (d)).

The $MnO_2$ and $Fe(OH)_{3A}$ concentrations matched the data points closely, however, the $Mn^{2+}$ and $Fe^{2+}$ concentrations in this simulation peaked at a higher concentration around 2 cm (Figure \@ref(fig:VCW-1) (e) to (h)). It was suspected that this might have resulted from poor resolution of the grid, however, this simulation had the same number and size of layers as Van Cappellen and Wang. Below 5 cm deep the simulation had a lower concentration than the field data, but the deepest 5 field data points for $Fe^{2+}$ were within the ranges of the simulations.

$H_2S$ had the greatest variation relative to its concentration, which increased with depth (Figure \@ref(fig:VCW-1) (i)). $FeS$ concentration had a consistently concave shape whereas the field data appeared to be concave around 3 cm deep, and the relatively high concentration of $FeS$ at 5 mmol L^-1^ in the oxic zone was not within the range of these simulations (Figure \@ref(fig:VCW-1) (j)). $FeS$ could be simulated close to the field data, but only by fixing the bottom boundary at 45 mmol L^-1^ and not because of its precipitation from $Fe^{2+}$ and $H_2S$.

```{r VCW-1, echo=FALSE, fig.cap="Concentration-depth profiles for major variables from 100 simulations, including the median result (black line), ranges of results (grey bands) and field data from Van Cappellen and Wang (1996) (squares).",fig.show='hold',fig.align='center', out.width = '80%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/VCW/VCW-01.png")
```

##### *Reaction rates* {.unnumbered}

The rates of aerobic respiration and denitrification using CANDI-AED were close to the simulated rates of Van Cappellen and Wang (1996), though in the deepest part of the sediment, the simulated aerobic rate was greater using this model (Figure \@ref(fig:VCW-2) (b)). As a result of the classic inhibition sequence, oxygen penetration to greater depths limited denitrification to a deeper depth layer (Figure \@ref(fig:VCW-2) (c)). This carried through to the peaks of all oxidation rates, which then occurred at slightly greater depths than in the Van Cappellen and Wang simulation. The maximum manganese reduction rate was smaller than in the Van Cappellen and Wang simulation, however, in both models the manganese reduction rate was very small relative to the other terminal oxidation processes (Figure \@ref(fig:VCW-2) (d)). The iron reduction rate was similar to the Van Cappellen and Wang rate, although slightly smaller, and could not be increased without leading to greater discrepancy in the dissolved iron concentration profiles (Figure \@ref(fig:VCW-2) (e)). The sulphate reduction profiles were very close and as with Van Cappellen and Wang, methanogenesis was completely inhibited (Figure \@ref(fig:VCW-2) (f)).

```{r VCW-2, echo=FALSE, fig.cap="Profiles of organic matter oxidation rate with depth. Plot (a) was the total organic matter reaction rate set by a fixed depth-oxidation rate, which Van Cappellen and Wang based on the measured data (squares). Plots (b) to (f) are these simulated rates (solid lines) compared with the simulated rates of Van Cappellen and Wang (dashed lines) and in the case of sulphate reduction, measured data (squares, (f))",fig.show='hold',fig.align='center', out.width = '80%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/VCW/R1.png")
```

### Case study: Aquaculture and sediment recovery

#### Introduction

This work was to calculate the potential effect of aquaculture on the sediment condition of fish cages propose for installation near Delma Island, UAE. The approach was similar to work done in 2015 on a proposed aquaculture site in the Abrolhos Islands, Western Australia (Paraska et al. 2015). The main objective of this part of the project was to calculate sediment chemical pore water concentrations and fluxes to the water column in response to variable deposition loads.

The primary input to the sediment model was reactive organic matter derived from aquaculture waste. This was calculated by models that tracked the transport and deposition of waste particles as a result of the hydrodynamics at the study site and different fish stocking scenarios. The organic matter deposition was then used as input to the sediment diagenesis model. The resulting waste particle distributions and sediment chemistry were then used to calculate the spatial effect of the cage operations on the sediment environment, as well as the greater chemistry and ecology of the water column.

The motivation behind the study was to ensure that potentially detrimental effects to the ecosystem were avoided. Fish cages may be positioned in one place for up to five years, after which the cage would be moved and the environment left to recover to its background condition (Figure \@ref(fig:UAEFish-1)).

```{r UAEFish-1, echo=FALSE, fig.cap="Schematic of the design of the aquaculture simulation. The simulation had a period of spin up, a period of aquaculture and then a recovery period.",fig.show='hold',fig.align='center', out.width = '80%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/UAEFish/UAE1-08.png")
```

Given enough energy in the water column, the concentrations of chemicals can return to background relatively quickly. However, the conditions in the sediment can take much longer to recover to background levels, because the transport of chemicals is much slower. In particular, there may be a high concentration of $H_2S$ and low concentration of $O_2$, which together would make the environment uninhabitable for benthic fauna. Further, there may be an ongoing flux of nitrogen and phosphorus from the sediment to the water column, which could fuel algal blooms.

#### Methods

##### *UAE study site chemical concentrations and fluxes* {.unnumbered}

Field data was measured at the study site in October 2017 and analysed in November 2017. CANDI-AED was tailored to as many local conditions as possible. Since the $O_2$ and $NO_3^-$ concentrations were high and the $NH_4^+$ concentration was zero, it was assumed that the water column was aerobic. Total nitrogen in the sediment was measured as approximately 300 mg kg^-1^, which was converted to approximately 30 mmol L^-1^. Inorganic phosphate and total phosphate in the bottom water were set to zero, since their concentrations were below the detection limit. Total phosphorus in the sediment was measured as approximately 500 mg kg^-1^, which was converted to approximately 20 mmol L^-1^.

Total organic carbon was measured as approximately 0.6% carbon g sediment^-1^. This was converted as approximately 600 mmol L^-1^, after dividing by molecular mass and multiplying by bulk density. The C:N:P ratio of background organic matter fluxes was set after calibrating the top 5 cm of sediment to field conditions. This was calibrated well to field conditions, however, it resulted in a ratio that is rich in N and P compared to refractory carbon that is typically present in sediments. It was decided that given this uncertainty, it would be best to calibrate to field conditions. This also resulted in a conservative approach to calculating background nutrient fluxes, in that they would not be underestimated. \$SO_4\^{2-} concentration was set as the same for the bottom water boundary and for the initial sediment porewater concentration, since this was a marine site.

Manganese and iron concentrations were below the detection limit in the water column field data, and so they were set to zero in the simulations. The reaction of iron and manganese with organic matter can affect the generation of $H_2S$, and therefore the ecological effect on the surface sediments. Therefore, an extra simulation was run with assumed low background iron and manganese influxes, to examine the effect, as outlined below. Manganese and iron concentrations in the sediment were measurable, and so to reconcile this difference between sediment and water concentrations, it was assumed that the measured sediment iron and manganese reflected crystallized, non-reactive minerals (known in diagenesis models as $Fe(OH)_{3B}$ and $MnO_{2B}$).

<center>

```{r param-UAE1, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
param_UAE1 <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet = "UAE" )
kable(param_UAE1 ,"html", escape = F, align = "c", caption = "Details of water chemistry specification for the sediment model",
      bootstrap_options = "hover") %>%
   kable_styling(param_Hidden, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "10em" , color = "black") %>%
 column_spec(2, width_min = "10em", color = "black") %>%
 column_spec(3, width_min = "10em", color = "black") %>%
 # column_spec(4, width_min = "5em") %>%
 row_spec(1:5, background = 'white', color = "black") %>%
 scroll_box(width = "45em", height = "20em",fixed_thead = FALSE)
```

</center>

##### *Aquaculture waste input parameters* {.unnumbered}

The C:N:P ratio of the aquaculture waste was supplied based on calculations based on a cage-model for fish biomass and waste properties. The rate constant for bacterial consumption of organic matter was set to 10 y^-1^. This is one of the most uncertain, yet crucial parameters in diagenesis modelling. The equivalent value used by Brigolin et al (2009) was 8 y^-1^. Considering that the water temperature at this study site may reach 40 °C, microbial metabolism, and therefore the rate constant, for all processes may approach a value twice as high as this. However, the value of 10 y^-1^ was used in order to take a conservative approach to environmental recovery time.

##### *Physical properties* {.unnumbered}

The simulated sediment was 20 cm deep, and had 100 layers of exponentially increasing size from the sediment-water interface (0.25 mm) to the deepest layer (4 mm). The rate of particle deposition to the sediment-water interface was set to 0.2 cm y-1. Sediment porosity was not measured at the study site and so an assumption of 0.6 (water/total volume) was assumed for the bulk of the sediment, increasing to 0.9 in the upper 5 cm.

##### *Simulation approach* {.unnumbered}

The overall approach was to assess the sediment over a "gradient" of different waste loading intensities. The setup was a matrix of forty simulations, with ten input fluxes tested, each for 4 aquaculture durations (table \@ref(tab:param-UAE2)).

Each simulation was for 23 years, with a 10 year period of background conditions, then the aquaculture period, then a recovery period. The duration of aquaculture was varied as either 1, 3 or 5 years. (For the larger fluxes, an additional period immediately following aquaculture was required, which had a fine time step.)

<center>

```{r param-UAE2, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
param_UAE2 <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet = "UAE2" )
kable(param_UAE1 ,"html", escape = F, align = "c", caption = "Details of water chemistry specification for the sediment model",
      bootstrap_options = "hover") %>%
   kable_styling(param_Hidden, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "10em" , color = "black") %>%
 column_spec(2, width_min = "15em", color = "black") %>%
 column_spec(3, width_min = "15em", color = "black") %>%
 # column_spec(4, width_min = "5em") %>%
 row_spec(1:5, background = 'white', color = "black") %>%
 scroll_box(width = "45em", height = "20em",fixed_thead = FALSE)
```

</center>

#### Results

##### *Fluxes* {.unnumbered}

Results of the steady state oxygen and nutrient fluxes being returned from the diagenesis model for the range of depositional rates (FOM) are shown below (figure \@ref(fig:UAEFish-2), table \@ref(tab:param-UAE4).

```{r UAEFish-2, echo=FALSE, fig.cap="Summary of the nutrient fluxes at the sediment water interface returned from the diagenesis model for a wide range of input dpeositional flux rates.",fig.show='hold',fig.align='center', out.width = '110%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/UAEFish/Results1.png")
```

\

<center>

```{r param-UAE4, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
param_UAE4 <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet = "UAE4" )
kable(param_UAE4 ,"html", escape = F, align = "c", caption = "Nutrient and oxygen flux results for several depositional categories.",
      bootstrap_options = "hover") %>%
   kable_styling(param_Hidden, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "7em" , color = "black", bold = T) %>%
 column_spec(2, width_min = "7em", color = "black") %>%
 column_spec(3, width_min = "7em", color = "black") %>%
 column_spec(4, width_min = "7em", color = "black") %>%
 column_spec(5, width_min = "7em", color = "black") %>%
 # column_spec(6, width_min = "7em", color = "black") %>%
 row_spec(1:5, background = 'white', color = "black") %>%
 scroll_box(width = "45em", height = "24em",fixed_thead = FALSE)
```

</center>

Recovery time periods are shown below in Figure \@ref(fig:UAEFish-3) UAEFish-3, and threshold depositional rates for recovery time categorisation is summarised in Table \@ref(tab:param-UAE6). The impact zone thresholds are summarised in Table \@ref(tab:param-UAE7).

\

```{r UAEFish-3, echo=FALSE, fig.cap="Predicted recovery time of the sediment for the 2cm oxygen and solfide concentration to return to within 85% of the pre-farmign value. Series depicted for aquaculture loading periods of 1-5 years.",fig.show='hold',fig.align='center', out.width = '90%'}
knitr::include_graphics("images/23-sediment_biogeochemistry/UAEFish/Results2.png")
```

<center>

```{r param-UAE6, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
param_UAE6 <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet = "UAE6" )
kable(param_UAE6 ,"html", escape = F, align = "c", caption = "Nutrient and oxygen flux results for several depositional categories.",
      bootstrap_options = "hover") %>%
   kable_styling(param_Hidden, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "12em" , color = "black", bold = T) %>%
 column_spec(2, width_min = "12em", color = "black") %>%
 column_spec(3, width_min = "12em", color = "black") %>%
 column_spec(4, width_min = "12em", color = "black") %>%
 column_spec(5, width_min = "12em", color = "black") %>%
 # column_spec(6, width_min = "7em", color = "black") %>%
 row_spec(1:5, background = 'white', color = "black") %>%
 scroll_box(width = "45em", height = "17em",fixed_thead = FALSE)
```

<br>

```{r param-UAE7, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(kableExtra)
library(readxl)
param_UAE7 <- read_xlsx("tables/23-sediment_biogeochemistry/UserManual.xlsx", sheet = "UAE7" )
kable(param_UAE7 ,"html", escape = F, align = "c", caption = "Depositional thresholds for four impact zone settings, summarised for 4 aquaculture operation periods.",
      bootstrap_options = "hover") %>%
   kable_styling(param_Hidden, bootstrap_options = "hover",
               full_width = F, position = "center",
               font_size = 12) %>%
 row_spec(0, background = "#14759e", bold = TRUE, color = "white") %>%
 column_spec(1, width_min = "7em" , color = "black", bold = T) %>%
 column_spec(2, width_min = "7em", color = "black") %>%
 column_spec(3, width_min = "7em", color = "black") %>%
 column_spec(4, width_min = "7em", color = "black") %>%
 # column_spec(5, width_min = "7em", color = "black") %>%
 # column_spec(6, width_min = "7em", color = "black") %>%
 row_spec(1:4, background = 'white', color = "black") %>%
 scroll_box(width = "45em", height = "12em",fixed_thead = FALSE)
```

</center>

\

### Case study: The Coorong Dynamics Model

CANDI-AED was used as part of the assessment of the Coorong Lagoon in South Australia. It was used as part of [a greater project](https://aquaticecodynamics.github.io/cdm-science/index.html), in conjunction with a fully 3D resolved Tuflow-FV model of the lagoon, that had full biogeochemistry modules and habitat assessments. In particular, the nutrient budgets and the effects of the hypersalinity were calculated. The [sediment modelling component of the project](https://aquaticecodynamics.github.io/cdm-science/sediment-biogeochemistry.html) in particular looked at the effect of hypersalinity on bioturbation and bioirrigation rates, and the resulting effects such as denitrification.

<!-- ```{r CoorongSVG-1, echo=FALSE, fig.cap="Test",fig.show='hold',fig.align='center', out.width = '90%'} -->
<!-- knitr::include_graphics("images/23-sediment_biogeochemistry/VCW/t_D_C_6panel_oxy.svg") -->
<!-- ``` -->

<!-- ```{r CoorongPNG-1, echo=FALSE, fig.cap="Test",fig.show='hold',fig.align='center', out.width = '100%'} -->
<!-- knitr::include_graphics("images/23-sediment_biogeochemistry/UAEFish/ThreePlots_popl_13_reso_500.png") -->
<!-- ``` -->

<!-- ```{r CoorongPNG-2, echo=FALSE, fig.cap="Test",fig.show='hold',fig.align='center', out.width = '100%'} -->
<!-- knitr::include_graphics("images/23-sediment_biogeochemistry/UAEFish/ThreePlots_popl_13_reso_1000.png") -->
<!-- ``` -->

<!-- ```{r CoorongPNG-3, echo=FALSE, fig.cap="Test",fig.show='hold',fig.align='center', out.width = '100%'} -->
<!-- knitr::include_graphics("images/23-sediment_biogeochemistry/UAEFish/ThreePlots_amm_11_reso_200.png") -->
<!-- ``` -->

## References

Berg, P., S. Rysgaard and B. Thamdrup (2003). Dynamic modeling of early diagenesis and nutrient cycling. A case study in an Arctic marine sediment. American Journal of Science 303(10): 905-955.

Berner, R. (1980). Early diagenesis: A theoretical approach. New Jersey, Princeton University Press.

Bethke, C., R. Sanford, M. Kirk, Q. Jin, and T. Flynn (2011), The thermodynamic ladder in geomicrobiology, American Journal of Science, 311(3): 183-210.

Boudreau, B. (1997). Diagenetic models and their implementation. Heidelberg, Springer.

Boudreau, B. and J. Westrich (1984). The dependence of bacterial sulfate reduction on sulfate concentration in marine sediments. Geochimica et Cosmochimica Acta 48(12): 2503-2516.

Boudreau, B. P. (1996). A method-of-lines code for carbon and nutrient diagenesis in aquatic sediments. Computers & Geosciences 22(5): 479-496.

Brigolin, D., R. Pastres, T. D. Nickell, C. J. Cromey, D. R. Aguilera and P. Regnier (2009). Modelling the impact of aquaculture on early diagenetic processes in sea loch sediments. Marine Ecology-Progress Series 388: 63-80.

Brown, P. N., G. D. Byrne, and A. C. Hindmarsh (1989), VODE: A variable-coefficient ODE solver, SIAM journal on scientific and statistical computing, 10(5): 1038-1051.

Dale, A. W., D. R. Aguilera, P. Regnier, H. Fossing, N. J. Knab and B. B. Jorgensen (2008a). Seasonal dynamics of the depth and rate of anaerobic oxidation of methane in Aarhus Bay (Denmark) sediments. Journal of Marine Research 66(1): 127-155.

Eldridge, P. and J. Morse (2008). Origins and temporal scales of hypoxia on the Louisiana shelf: Importance of benthic and sub-pycnocline water metabolism. Marine Chemistry 108(3-4): 159-171.

Fossing, H., P. Berg, B. Thamdrup, S. Rysgaard, H. Sorensen and K. Nielsen (2004). A model set-up for an oxygen and nutrient flux model for Aarhus Bay (Denmark).

Gaillard, J. F. and C. Rabouille (1992). Using Monod kinetics in geochemical models of organic carbon mineralisation in deep-sea surficial sediments. Deep-Sea Food Chains and the Global Carbon Cycle. G. T. Rowe and V. Pariente. 360: 309-324.

Jin, Q. S. and C. M. Bethke (2005). Predicting the rate of microbial respiration in geochemical environments. Geochimica et Cosmochimica Acta 69(5): 1133-1143.

Luff, R., K. Wallmann, S. Grandel and M. Schluter (2000). Numerical modeling of benthic processes in the deep Arabian Sea. Deep-sea research. Part 2, Topical studies in oceanography 47(14): 3039-3072.

Reed, D. C., C. P. Slomp and B. G. Gustafsson (2011b). Sedimentary phosphorus dynamics and the evolution of bottom-water hypoxia: A coupled benthic-pelagic model of a coastal system. Limnology and Oceanography 56(3): 1075-1092.

Soetaert, K., P. M. J. Herman and J. J. Middelburg (1996a). A model of early diagenetic processes from the shelf to abyssal depths. Geochimica et Cosmochimica Acta 60(6): 1019-1040.

Tufano, K. J., S. G. Benner, K. U. Mayer, M. A. Marcus, P. S. Nico, and S. Fendorf (2009), Aggregate-Scale Heterogeneity in Iron (Hydr)oxide Reductive Transformations All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher, Vadose Zone J., 8(4): 1004-1012.

Van Cappellen, P. and Y. F. Wang (1996). Cycling of iron and manganese in surface sediments: A general theory for the coupled transport and reaction of carbon, oxygen, nitrogen, sulfur, iron, and manganese. American Journal of Science 296(3): 197-243.

<!-- install.packages("kableExtra", lib="C:/Program Files/R/R-4.2.0/library") -->

<!-- install.packages("knitr", lib="C:/Program Files/R/R-4.2.0/library") -->

<!-- install.packages("readxl", lib="C:/Program Files/R/R-4.2.0/library") -->

<!-- library(rmarkdown, lib.loc = "C:/Program Files/R/R-4.2.0/library") -->

<!-- install.packages("downlit") -->

<!-- library(knitr) -->

<!-- library(kableExtra) -->

<!-- library(readxl) -->

<!-- library(downlit) -->

<!-- <!-- @ref(fig:23-pic2) --> -->

<!-- install.packages("hablar") -->

<!-- install.packages("xaringanextra") -->

<!-- install.packages("leaflet") -->

<!-- install.packages("gt") -->

<!-- install.packages("seacarb") -->

<!-- install.packages("plotly") -->

<!-- install.packages("tidyverse") -->