running-the-model.qmd

---
title: Process rate estimation
execute: 
  eval: false
  echo: true
---

## State function set

The process rate estimator includes three-state functions, describing the change in N~2~O concentrations over time.
The change in N~2~O concentrations in each depth increment over time $\frac{\Delta}{\Delta t}[\ce{N2O}]$ depends on the flux of N~2~O entering the depth increment from the top or the bottom through diffusion ($\text F_\text{t,i}$ and $\text F_\text {b,i}$, respectively), the flux of N~2~O leaving the depth increment through diffusion ($\text F_\text {out}$), the rate of N~2~O produced through nitrification (`N2Onit`), the rate of N~2~O produced through denitrification (`N2Oden`), and the rate of N~2~O reduced to N~2~ (`N2Ored`).

$$\frac{\Delta[\ce{N2O}]}{\Delta t} = \text F_{\text{t,i}} + \text F_{\text{b,i}}  + \ce{N2O}_{\text{nit.}} + \ce{N2O}_{\text{den.}} + \ce{N2O}_{\text{red.}}$${#eq-stateN2O}

For the equation of the change of site preference over time, we'll use the substitutions $A_1$, $A_2$ and $A_3$, as the expression is rather long.

$$\begin{gathered}
\frac{\Delta \text{SP}}{\Delta t} = \frac{A_1 + A_2 + A_3}{[\ce{N2O}]} \\ \\
A_1 = \text F_{\text{t,i}}(\text{SP}_{\text{t,i}} - \eta\: \text{SP}_\text{diff.} - \text{SP}) + \text F_{\text{b,i}}(\text{SP}_{\text{b,i}} - \eta\: \text{SP}_\text{diff.} - \text{SP}) \\
A_2 = \ce{N2O}_{\text{nit.}}(\text{SP}_\text{nit.} - \text{SP}) + \ce{N2O}_{\text{den.}}(\text{SP}_\text{den.} - \text{SP}) \\
A_3 = - \eta\: \text{SP}_\text{diff.} \text F_{\text{out}} - \eta\text{SP}_\text{red.} \ce{N2O}_{\text{red.}}
\end{gathered}$${#eq-stateSP}
Similar to $\text{SP}$, the change in $\delta^{18}\text O$ over time in each time point and depth increment can be described as:

$$\begin{gathered}
\frac{\Delta\:\delta^{18}\text{O}}{\Delta t} = \frac{B_1 + B_2 + B_3}{[\ce{N2O}]} \\ \\
B_1 = \text F_{\text{t,i}}(\delta^{18}\text{O}_{\text{t,i}} - \eta\: ^{18}\text{O}_\text{diff.} - \delta^{18}\text{O}) + \text F_{\text{b,i}}(\delta^{18}\text{O}_{\text{b,i}} - \eta\: ^{18}\text{O}_\text{diff.} - \delta^{18}\text{O}) \\
B_2 = \ce{N2O}_{\text{nit.}}(\delta^{18}\text{O}_\text{nit.} - \delta^{18}\text{O}) + \ce{N2O}_{\text{den.}}(\delta^{18}\text{O}_\text{den.} - \delta^{18}\text{O}) \\
B_3 = - \eta\: ^{18}\text{O}_\text{diff.} \text F_{\text{out}} - \eta\ce{^{18}O_\text{red.}} \ce{N2O}_{\text{red.}}
\end{gathered}$${#eq-stated18O}

## How to run the process rate estimator in R

To run the process rate estimator in R, you'll first have to install the package `PRE` as described [in the introduction](https://damian-oswald.github.io/process-rate-estimator/introduction.html#installation-of-the-r-package-pre).
Once installed, add the package to the search path.

```{r}
library(PRE)
```

### Prepare the data

In a first step, we can calculate the N~2~O-N from the original measurements.

```{r}
data <- getN2ON(data = PRE::measurements, parameters = PRE::getParameters())
```

Next, the data are interpolated. This process uses kernel regression to interpolate between the measured values of N~2~O-N, SP and d18O.

```{r}
data <- getMissing(data = data, hyperparameters = PRE::hyperparameters)
```

With all of these steps done, we can calculate the Fluxes according to equations 1, 2, 3.

```{r}
data <- calculateFluxes(data = data, parameters = PRE::getParameters())
```

Now, everything is ready to run the PRE.

### Running the process rate estimator

```{r}
# run the solver once
x <- PRE(data = data, column = 1, depth = 7.5, date = "2016-01-01", n = 1000)

# print out information
print(x)

# plot
plot(x)

# pairs panel
pairs(x)
```

Run the process rate estimator over the entire time span.

```{r}
x = longPRE(data, column = 1, depth = 7.5, n = 15)
```

Run the estimation of process rates for the entire project, i.e. applying the `longPRE` function to all combinations of depths and columns.
For very robust results, `n = 1000` is used in the corresponding R script.

```zsh
Rscript scripts/run-process-rate-estimator/run-PRE.R
```

The results will be written as `estimated-process-rates.csv` into the `output` directory.