---
title: "Getting Started with okBATHTUB"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting Started with okBATHTUB}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 7,
  fig.height = 4.5,
  dpi       = 120
)
```

## Overview

`okBATHTUB` is an R implementation of steady-state empirical reservoir
eutrophication modelling. It supports three retention model families:

1. **Walker (1985, 1996) BATHTUB Model 1** — the default — a
   second-order available-phosphorus sedimentation model calibrated to
   U.S. Army Corps of Engineers reservoir data.
2. **Vollenweider (1976) / Larsen-Mercier (1976)** — a parsimonious
   first-order hydraulic-residence retention model. Equivalent to
   Walker BATHTUB Model 5 ("Northern Lakes"), which Walker (1996)
   explicitly flags as *not* calibrated to Corps of Engineers reservoir
   data.
3. **Oklahoma** — Walker Model 1 retention paired with Oklahoma-
   specific chlorophyll-a and Secchi regression coefficients calibrated
   from publicly available state lake monitoring data.

The package predicts in-lake total phosphorus, total nitrogen,
chlorophyll-a, and Secchi depth, and computes Carlson (1977) Trophic
State Indices. It is designed to complement watershed loading models
such as SWAT or HAWQS in a two-model nutrient management workflow.

```{r load}
library(okBATHTUB)
```

## The core pipeline

The standard single-segment workflow chains five functions with the
base pipe operator:

```
ok_load() |> ok_hydraulics() |> ok_retention() |> ok_inlake() |> ok_tsi()
```

### Step 1 — Load inputs

`ok_load()` accepts tributary inflow volume and flow-weighted nutrient
concentrations:

```{r load_step}
result <- ok_load(
  inflow_m3yr   = 45e6,    # tributary inflow (m^3/yr)
  tp_inflow_ugl = 120,     # flow-weighted mean TP (ug/L)
  tn_inflow_ugl = 1800     # flow-weighted mean TN (ug/L)
)
print(result)
```

### Step 2 — Hydraulics

`ok_hydraulics()` adds reservoir morphometry and computes residence
time and areal water load:

```{r hydraulics_step}
result <- result |>
  ok_hydraulics(
    surface_area_ha = 890,
    mean_depth_m    = 4.2
  )

cat("Hydraulic residence time:", round(result$data$hydraulic_residence_time_yr, 2), "yr\n")
cat("Areal water load (qs):   ", round(result$data$areal_water_load_myr, 2), "m/yr\n")
```

### Step 3 — Nutrient retention

`ok_retention()` estimates the fraction of incoming TP and TN retained
in the reservoir. The default uses Walker BATHTUB Model 1
(second-order available-P sedimentation):

$$P = \frac{-1 + \sqrt{1 + 4 K A_1 P_i \tau}}{2 K A_1 \tau}$$

where $A_1 = 0.17 \cdot Q_s/(Q_s + 13.3)$ and
$Q_s = \max(Z/\tau, 4)$ m/yr. The retention coefficient stored on the
object is back-calculated from this solution so that the downstream
mass balance $C_{lake} = C_{in}(1 - R)$ exactly reproduces the
Model 1 result.

```{r retention_step}
result <- result |> ok_retention()

cat("TP retention coefficient:", round(result$data$tp_retention_coeff, 3),
    "(", result$data$tp_retention_form, ")\n")
cat("TN retention coefficient:", round(result$data$tn_retention_coeff, 3), "\n")
```

To use the simpler Vollenweider / Larsen-Mercier form instead, pass
`coefficients = "vollenweider"` to `ok_load()`.

### Step 4 — In-lake predictions

`ok_inlake()` applies the mass balance and the chlorophyll-a / Secchi
regressions:

```{r inlake_step}
result <- result |> ok_inlake()

cat("In-lake TP:    ", round(result$data$tp_inlake_ugl, 1), "ug/L\n")
cat("Chlorophyll-a: ", round(result$data$chla_ugl, 2),    "ug/L\n")
cat("Secchi depth:  ", round(result$data$secchi_m, 2),    "m\n")
```

### Step 5 — Trophic State Index

`ok_tsi()` computes Carlson (1977) TSI from all predicted variables:

```{r tsi_step}
result <- result |> ok_tsi()
summary(result)
```

### Full pipeline in one call

```{r full_pipeline}
result <- ok_load(
    inflow_m3yr   = 45e6,
    tp_inflow_ugl = 120,
    tn_inflow_ugl = 1800
  ) |>
  ok_hydraulics(surface_area_ha = 890, mean_depth_m = 4.2) |>
  ok_retention() |>
  ok_inlake()    |>
  ok_tsi()

summary(result)
```

## Using the bundled reservoir lookup

The `ok_reservoirs` dataset contains morphometry for major Oklahoma
reservoirs compiled from publicly available USACE, BOR, NID, and OWRB
sources. Use `ok_reservoir()` for a quick lookup:

```{r reservoir_lookup}
res <- ok_reservoir("Arcadia Lake", exact = TRUE)
res[, c("lake_name", "surface_area_ha", "mean_depth_m", "eco_l3_name",
        "data_quality")]
```

```{r reservoir_pipeline}
result <- ok_load(
    inflow_m3yr   = 45e6,
    tp_inflow_ugl = 120,
    tn_inflow_ugl = 1800
  ) |>
  ok_hydraulics(
    surface_area_ha = res$surface_area_ha[1],
    mean_depth_m    = res$mean_depth_m[1]
  ) |>
  ok_retention() |>
  ok_inlake()    |>
  ok_tsi()

cat("Trophic state:", result$data$trophic_state, "\n")
cat("Mean TSI:     ", round(result$data$tsi_mean, 1), "\n")
```

For a dataset overview by ecoregion:

```{r reservoir_summary}
ok_reservoir_summary()
```

## Multiple tributaries

When more than one tributary contributes to the reservoir, use
`ok_load_multi()` to compute flow-weighted mean concentrations
automatically:

```{r multi_trib}
tributaries <- data.frame(
  inflow_m3yr   = c(30e6, 15e6),
  tp_inflow_ugl = c(100,  160),
  tn_inflow_ugl = c(1500, 2400)
)

result_multi <- ok_load_multi(tributaries) |>
  ok_hydraulics(surface_area_ha = 890, mean_depth_m = 4.2) |>
  ok_retention() |>
  ok_inlake()    |>
  ok_tsi()

cat("FW mean TP inflow:", round(result_multi$data$tp_inflow_ugl, 1), "ug/L\n")
```

## Computing TSI from observed data

`ok_tsi_observed()` computes Carlson TSI directly from monitoring data
without running the full pipeline — useful for comparing observed
trophic state against modelled predictions:

```{r tsi_observed}
obs <- ok_tsi_observed(
  tp_ugl    = 85,
  chla_ugl  = 22,
  secchi_m  = 0.8
)

cat("TSI(TP):    ", round(obs$tsi_tp,     1), "\n")
cat("TSI(Chl-a): ", round(obs$tsi_chla,   1), "\n")
cat("TSI(Secchi):", round(obs$tsi_secchi, 1), "\n")
cat("Mean TSI:   ", round(obs$tsi_mean,   1), "\n")
cat("Class:      ", obs$trophic_state,         "\n")
```

## Comparing coefficient sets

A useful sanity check is to run the same reservoir under different
coefficient assumptions and see how predictions vary:

```{r coeff_comparison}
run_pipeline <- function(coef, eco = NULL) {
  ok_load(
      inflow_m3yr   = 45e6,
      tp_inflow_ugl = 120,
      tn_inflow_ugl = 1800,
      coefficients  = coef,
      ecoregion     = eco
    ) |>
    ok_hydraulics(surface_area_ha = 890, mean_depth_m = 4.2) |>
    ok_retention() |>
    ok_inlake()    |>
    ok_tsi()
}

r_walker  <- run_pipeline("walker")
r_voll    <- run_pipeline("vollenweider")
r_okla_ct <- run_pipeline("oklahoma", "Cross Timbers")

comparison <- data.frame(
  scenario   = c("Walker Model 1",
                 "Vollenweider/Larsen-Mercier",
                 "Oklahoma (Cross Timbers)"),
  tp_inlake  = round(c(r_walker$data$tp_inlake_ugl,
                       r_voll$data$tp_inlake_ugl,
                       r_okla_ct$data$tp_inlake_ugl), 1),
  chla       = round(c(r_walker$data$chla_ugl,
                       r_voll$data$chla_ugl,
                       r_okla_ct$data$chla_ugl), 2),
  tsi_mean   = round(c(r_walker$data$tsi_mean,
                       r_voll$data$tsi_mean,
                       r_okla_ct$data$tsi_mean), 1)
)
print(comparison, row.names = FALSE)
```

A few observations to expect:

- Walker Model 1 typically predicts **higher retention** (lower in-lake
  TP) than Vollenweider for the same residence time, because the
  second-order kinetics scale settling with concentration.
- Oklahoma Cross Timbers Chl-a coefficients produce a **shallower
  log-log slope** than Walker's national regression: for high TP,
  Oklahoma predicts *less* algal response per unit phosphorus,
  consistent with the inorganic-turbidity-dominated nature of
  many Cross Timbers reservoirs.

## References

- Carlson, R.E. (1977). A trophic state index for lakes.
  *Limnology and Oceanography*, 22(2), 361-369.
- Larsen, D.P. & Mercier, H.T. (1976). Phosphorus retention capacity
  of lakes. *Journal of the Fisheries Research Board of Canada*, 33(8),
  1742-1750.
- Vollenweider, R.A. (1976). Advances in defining critical loading
  levels for phosphorus in lake eutrophication. *Memorie dell'Istituto
  Italiano di Idrobiologia*, 33, 53-83.
- Walker, W.W. (1985). Empirical methods for predicting eutrophication
  in impoundments; Report 3, Phase III: Model refinements. Technical
  Report E-81-9, U.S. Army Engineer Waterways Experiment Station.
- Walker, W.W. (1996). *Simplified Procedures for Eutrophication
  Assessment and Prediction: User Manual*. Instruction Report W-96-2,
  U.S. Army Engineer Waterways Experiment Station.