module2.Rmd

# Module 2: Biogeochemistry {#M2_1 -}

<h2>Bacterially Mediated Organic Matter Oxidation Rates</h2>

### Objectives {#M2_2 -}
Become familiar with basic spread sheet commands in Excel by building a simple biogeochemical model from chemical and mathematical equations. Become familiar with the program R Studio and compare it to Excel.

### Before you start {#M2_before_start -}

- Listen to lecture on ordinary differential equations (ODEs)
- Be familiar with the exponential growth of Quokkas example in Excel
- Review slides on sediment

### Modelling concepts {#M2_3 -}
This model has variables of concentration and rate, no spatial dimensions, and it has a dimension of time. Last week’s model had variables of volume and concentration (salinity), and was also resolved in time. In coming weeks the models will be resolved in zero, one, two or three dimensions of space, and only some will be resolved in time.

Remember that a chemical reaction can be described using a simple differential equation. For the concetration of OM, we can write it as:
<center>
<br>
\begin{equation}
\frac{dOM}{dt} = -R_{OM} =  -k_{OM}OM
(\#eq:module2-1a)
\end{equation}
</center>
<br>
Here, $k_{OM}$ is the rate of the organic matter decay reaction, in units of $/year$, and $OM$ is a concentration of organic matter in units of $mM$.

In this activity going down the rows of the Excel spread sheet means going forward in time, as the differential equations divide the time axis into hundreds of tiny little chunks (i.e., each row is a day in this weeks setup). The concentration at any point in time may depend on concentration at the previous point in time. We use the notation concentration @ ^t^ depends on concentration @ ^t-1^. We write this by "discretizing" the equation, as:
<br>
\begin{equation}
OM^{t} = OM^{t-1} -  (k_{OM} OM^{t-1})\Delta t
(\#eq:module2-1b)
\end{equation}
</center>
<br>
This states that the current concentration is equal to the previous days concentration, subtract any reaction that occured over the past day.

### Organic matter oxidation {#M2_4 -}
Organic matter (dead algae, leaves etc.) is oxidised by different bacteria that use a range of oxidants. The bacteria simultaneously consume organic matter and an oxidant, in a similar way to how humans consume carbohydrate and oxygen. Scientists have observed that there is a general sequence of oxidants used, based on the energy available for each process. When all of the oxidants are consumed, methanogenesis takes place, where bacteria use the organic matter alone without any oxidant. This is similar to the lactic acid metabolic system that the human body uses when oxygen is in short supply, such as when you go for a run. Remember that methane is a reduced byproduct, not an oxidant. Other reduced by-products include N~2~ gas, Mn^2+^ and Fe^2+^ metal ions, and the S^2-^ anion.

We can construct a model to determine the rates at which the oxidation takes place, and use the rates to determine other processes such as oxygen consumption, denitrification or phosphate adsorption.

The rate is dependent on the concentration of the organic matter, $\color{#E36C0A}{OM}$, and a rate constant, $\color{#00B0F0}{k_{OM}}$. The rate is not dependent on the concentration of the oxidant, until the oxidant concentration becomes very low: this process is included as a <span style="color: #76923C;">limitation</span> term (With more O~2~ you do not become super; with less O~2~ you get sleepy). A <span style="color: #76923C;">Monod</span> equation is used to limit the rate.

<center>
<br>
\begin{equation}
R_{O_{2}} = \color{#00B0F0}{k_{OM}}\color{#E36C0A}{OM}\color{#76923C}{\frac{O_{2}}{O_{2}+K_{O_{2}}}}
(\#eq:module2-1)
\end{equation}
</center>
<br>

```{r M2Plot1, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="Move the slider to see the effect of changing the $K_{O_2}$ value in the monod equaiton.", out.width = '100%'}
library(plotly)
library(dplyr)

monod <- data.frame(
  O2 = double(),
  KO2 = double(),
  O2KO2 = double()
)
for(i in seq(5, 100, by = 5)){
  monod_temp <- data.frame(
    O2 = seq(0, 50, by = 5),
    KO2 = i
  )
  monod_temp$O2KO2 <- (monod_temp$O2)/(monod_temp$O2+monod_temp$KO2)  
  monod <- bind_rows(monod,monod_temp)
}


fig <- plot_ly()
fig <- fig %>% add_trace(
  monod,
  name="Monod",
  y = monod$O2KO2,
  x = monod$O2,
  frame = monod$KO2,
  line = list(color='rgb(118, 146, 60)'),
  hovertemplate = '%{y:.1f}<extra></extra>',
  type = 'scatter',
  mode = 'lines')
fig <- fig %>%
  layout(
    xaxis = list(title='O<sub>2</sub> Concentration'), 
    yaxis = list(title='O<sub>2</sub>/(O<sub>2</sub>+K<sub>O<sub>2</sub></sub>)')
  ) %>% 
  animation_opts(
    750, easing = "linear", redraw = FALSE
  ) %>% 
  animation_slider(
    currentvalue = list(prefix = "K<sub>O<sub>2</sub></sub>: ")
  ) %>% 
  config(displayModeBar = FALSE,
         displaylogo = FALSE)
fig  
```

The generally observed sequence of oxidants is modelled with an <span style="color: #FF0000;">inhibition</span> process, so that a higher energy yielding oxidant, such as oxygen, inhibits a lower energy oxidant such as iron.

<center>
<br>
\begin{align*}
R_{NO_{3}} &= \color{#00B0F0}{k_{OM}}\color{#E36C0A}{OM}\color{#76923C}{\frac{NO_{3}}{NO_{3}+K_{NO_{3}}}}\color{#FF0000}{\frac{K_{O_{2}}}{O_{2}+K_{O_{2}}}}
(\#eq:module2-2)\\
R_{Mn} &= \color{#00B0F0}{k_{OM}}\color{#E36C0A}{OM}\color{#76923C}{\frac{Mn}{Mn+K_{Mn}}}\color{#FF0000}{\frac{K_{NO_{3}}}{NO_{3}+K_{NO_{3}}}}\color{#FF0000}{\frac{K_{O_{2}}}{O_{2}+K_{O_{2}}}}
(\#eq:module2-3)\\
R_{Fe} &= \color{#00B0F0}{k_{OM}}\color{#E36C0A}{OM}\color{#76923C}{\frac{Fe}{Fe+K_{Fe}}}\color{#FF0000}{\frac{K_{Mn}}{Mn+K_{Mn}}}\color{#FF0000}{\frac{K_{NO_{3}}}{NO_{3}+K_{NO_{3}}}}\color{#FF0000}{\frac{K_{O_{2}}}{O_{2}+K_{O_{2}}}}
(\#eq:module2-4) \\
R_{SO_{4}} &= \color{#00B0F0}{k_{OM}}\color{#E36C0A}{OM}\color{#76923C}{\frac{SO_{4}^{2-}}{SO_{4}^{2-}+K_{SO_{4}}}}\color{#FF0000}{\frac{K_{Fe}}{Fe+K_{Fe}}}\color{#FF0000}{\frac{K_{Mn}}{Mn+K_{Mn}}}\color{#FF0000}{\frac{K_{NO_{3}}}{NO_{3}+K_{NO_{3}}}}\color{#FF0000}{\frac{K_{O_{2}}}{O_{2}+K_{O_{2}}}}
(\#eq:module2-5)\\
R_{Meth} &= \color{#00B0F0}{k_{OM}}\color{#E36C0A}{OM}\color{#FF0000}{\frac{K_{SO_{4}}}{SO_{4}^{2-}+K_{SO_{4}}}}\color{#FF0000}{\frac{K_{Fe}}{Fe+K_{Fe}}}\color{#FF0000}{\frac{K_{Mn}}{Mn+K_{Mn}}}\color{#FF0000}{\frac{K_{NO_{3}}}{NO_{3}+K_{NO_{3}}}}\color{#FF0000}{\frac{K_{O_{2}}}{O_{2}+K_{O_{2}}}}
(\#eq:module2-6)\\
\end{align*}
</center>
<br>

Note that the total amount of $OM$ decay, $R_{OM}$, is the sum of the all the individual pathways listed above.

In this lab, we simulate a blob of aquatic sediment. It is a closed system, so that the initial concentrations react until they reach equilibrium. This system is also perfectly mixed, with an unchanging biomass of bacteria, and so there is no spatial resolution. This allows us to focus on the chemical reactions alone.

### Module resources {#M2_5 -}

Download the Excel spreadsheet and R scripts for this module by clicking the download button in the tool bar <i class="fa fa-download" aria-hidden="true"></i>.

### Exercises {#M2_6 -}

#### Building the model in Excel {#M2_7 -}

1) Looking at the sequence of equations above: what is the difference between the limitation and the inhibition term? Have a good think about this. If concentration were amazingly bigger or smaller than the $K_{O_{x}}$^[Where $K_{O_{x}}$ is the general term for the half saturation constant of any oxidant, for example $K_{O_{2}}$, $K_{SO_{4}}$ etc.], what would happen to the rates of the higher and lower energy processes? For example, with a high oxygen concentration, what would happen to the manganese reduction rate? With a low oxygen and nitrate concentration, and a high manganese concentration, what would happen to the manganese oxidation rate?


2) Complete the Excel spreadsheet for iron reduction, sulfate reduction and methanogenesis. Oxygen, nitrate and manganese are done already, so copy the method and use the equations at the tops of the columns. It is tedious, fiddly and there is a good chance that you will introduce bugs that are hard to find. This is the nature of Excel.

```{r M2Vid1, echo=FALSE, message=FALSE, warning=FALSE}
htmltools::tags$video(
  width='100%', 
  controls=NA,
  htmltools::tags$source(
    src='images/module2/video1.mp4',
    type="video/mp4"
  )
  )
```

```{block2, hintM2-1, type='rmdnote2'}
The rate is dependent on the concentration in the previous time step – this is called an *explicit* solution (as opposed to *implicit*).
```

3) Complete the line plot (**1**) of relative oxidant concentration for iron and sulphate. Right click on the plot and use ‘Add data’, make the x axis time and the y axis relative concentration. It is easy to visualise if you keep the colours consistent, and choose the colours that you like. Adjust the values for the initial condition and the half-saturation constants and see how the lines move on your plot. Add more time steps underneath if you want to see what happens over the long term.

4) Complete the area plot (**2**) of reaction rates, then play around with the numbers to test the sensitivities of the outputs to the parameter inputs

```{r M2Plot2, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="Move the slider to experiment with changing the $K_{O_2}$ value.", out.width = '100%'}
library(plotly)
library(dplyr)
monod <- readRDS('data/monodData.rds')
fig <- plot_ly()
fig <- fig %>% add_trace(
  monod,
  name="OM",
  x = monod$time,
  y = monod$OMprop,
  frame = monod$KO2,
  line = list(color='rgb(227, 108, 10)'),
  hovertemplate = 'OM: %{y:.2f}<extra></extra>',
  type = 'scatter',
  mode = 'lines')
fig <- fig %>% add_trace(
  monod,
  name="O<sub>2</sub>",
  x = monod$time,
  y = monod$O2prop,
  frame = monod$KO2,
  line = list(color='rgb(118, 146, 60)'),
  hovertemplate = 'O<sub>2</sub>: %{y:.2f}<extra></extra>',
  type = 'scatter',
  mode = 'lines')
fig <- fig %>%
  layout(
    xaxis = list(title='Time',hoverformat = "Time"), 
    yaxis = list(title='Relative oxidant concentration')
  ) %>% 
  animation_opts(
    750, easing = "linear", redraw = FALSE
  ) %>% 
  animation_slider(
    currentvalue = list(prefix = "K<sub>O<sub>2</sub></sub>: ")
  ) %>% 
  config(displayModeBar = FALSE,
         displaylogo = FALSE)
fig  
```

####  Building the model in R {#M2_8 -}

Now the Excel spreadsheet is working, the next step is to repeat the process in R.

5) Open the R project and the three scripts: there is one in which you set the input parameters and initial conditions; one which has the equations; and one which is used for plotting the outputs. Read through them and become familiar with what is in each.

6) Without making any edits, run each script in sequence, starting with Inputs, then Equations, then Plotting. The plotting script should make a plot appear in the bottom right corner.

```{block2, hintM2-2, type='rmdtip2'}
You can press `Cmd+Enter` on MacOS, or `Ctrl+Enter` on Windows instead of clicking the 'Run' button.
```

7) You'll  have noticed some parts of these scripts are incomplete (i.e. commented out with `#`s). Complete the model by entering the inputs, equations and plotting scripts for iron, sulfate and methanogenesis reactions, just as you did in the Excel spreadsheet. Remember that you just need a rate of methanogenesis: methane is not an oxidant like sulfate.

::::: {.rmdtip2}

Step 1: Enter the inputs for Iron. Don't forget to 'run' these!

```{r M2Vid2, echo=FALSE, message=FALSE, warning=FALSE}
htmltools::tags$video(
  width='100%', 
  controls=NA,
  htmltools::tags$source(
    src='images/module2/video2.mp4',
    type="video/mp4"
  )
  )
```
<br>
Step 2: Write the equation. It should match Equation \@ref(eq:module2-4).

```{r M2Vid3, echo=FALSE, message=FALSE, warning=FALSE}
htmltools::tags$video(
  width='100%', 
  controls=NA,
  htmltools::tags$source(
    src='images/module2/video3.mp4',
    type="video/mp4"
  )
  )
```
<br>
Step 3: Graph it and repeat previous steps for the other oxidants.

```{r M2Vid4, echo=FALSE, message=FALSE, warning=FALSE}
htmltools::tags$video(
  width='100%', 
  controls=NA,
  htmltools::tags$source(
    src='images/module2/video4.mp4',
    type="video/mp4"
  )
  )
```
:::::

8) In the Plotting script, graph the concentrations, and relative concentrations, of each oxidant using `ggplot()`

Now we have a working process-based model. This was the hard and boring part, but now we can be creative. So far we have established an initial condition and let it run to completion, but next we can model a hypothetical environmental situation. For example, imagine there were an oxygen fluctuation, perhaps from algal photosynthesis occurring in the day but not at night in a lake, or from the tidal movement of hypoxic water.

9) In either Excel or R, create a new environmental situation that involves organic matter oxidation. Consider the assumptions you are making and how important they are. Use the model to test the effect of this situation on concentrations and rates. For example, what would happen if there were regular organic matter deposition as a boundary condition? Or maybe what if there were a constant oxygen concentration? What conditions might be needed to minimise methanogenesis? Be creative! Write a brief description of your new situation and create one plot (**4**) to explain your result.

### Submission {#M2_9 -}
:::: {.redbox2}

**Submit properly formatted graphs and tables of the following sections of the lab:**

- Graph of relative and absolute concentration (y) against time (x) scatter plot, for all oxidants, made in Excel
- Graph of rate (y) against time (x) area plot, for all oxidants
- Graph of relative and absolute concentration (y) against time (x) scatter plot, for all oxidants, made in R
- Any figure of a new environmental situation you create, and a brief explanation of less than 100 words of how you used the model creatively.

These are to be uploaded as per the formatting specified in the [Style Guide](#StyleGuide). Marks will be deducted for incorrect formatting.
<!--
These are in a word doc or PDF format. No screenshots of figures from Excel/Excel spreadsheets to be uploaded

**General professional formatting guidelines:**

- All figures are to have adequate captions explaining them
- For graphs, figure captions go below the plot
- For tables, the caption goes above the table
- Make sure figures and their text size is readable

**Excel hints:**

- When there is a caption for a plot, you remove the title
- Remove the plot border and gridlines
- Make sure both axes have visible lines and tick marks
- Units need to be noted properly with the axis label - 'Temperature (°C)'
- Round numbers to be reasonable
-->
::::