Getting Started with okBATHTUB

Overview

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

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.
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.
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.

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:

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)
#> -- okBATHTUB Result --
#>   Pipeline step : load
#>   Segment       : main
#>   Coefficients  : walker
#> 
#>   inflow_m3yr    : 4.5e+07
#>   tp_inflow_ugl  : 120
#>   tp_load_kgyr   : 5400
#>   tn_inflow_ugl  : 1800
#>   tn_load_kgyr   : 8.1e+04

Step 2 — Hydraulics

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

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")
#> Hydraulic residence time: 0.83 yr
cat("Areal water load (qs):   ", round(result$data$areal_water_load_myr, 2), "m/yr\n")
#> Areal water load (qs):    5.06 m/yr

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.

result <- result |> ok_retention()

cat("TP retention coefficient:", round(result$data$tp_retention_coeff, 3),
    "(", result$data$tp_retention_form, ")\n")
#> TP retention coefficient: 0.632 ( walker_model1 )
cat("TN retention coefficient:", round(result$data$tn_retention_coeff, 3), "\n")
#> TN retention coefficient: 0.553

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:

result <- result |> ok_inlake()

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

Step 5 — Trophic State Index

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

result <- result |> ok_tsi()
summary(result)
#> ========================================
#>   okBATHTUB Water Quality Summary
#> ========================================
#> 
#>   Segment      : main
#>   Coefficients : walker
#>   Pipeline     : tsi
#> 
#>   -- Hydraulics --
#>   Inflow           : 4.500e+07 m3/yr
#>   Surface area     : 890.0 ha
#>   Mean depth       : 4.20 m
#>   Residence time   : 0.831 yr
#>   Areal water load : 5.06 m/yr
#> 
#>   -- Nutrient Retention --
#>   TP retention     : 0.632  (walker_model1)
#>   TN retention     : 0.553  (walker_model1)
#> 
#>   -- In-Lake Predictions --
#>   TP               : 44.2 ug/L
#>   TN               : 803.8 ug/L
#>   Chlorophyll-a    : 17.68 ug/L
#>   Secchi depth     : 1.06 m
#> 
#>   -- Carlson Trophic State Index --
#>   TSI(TP)          : 58.8
#>   TSI(Chl-a)       : 58.8
#>   TSI(Secchi)      : 59.1
#>   TSI(mean)        : 58.9  (n = 3 components)
#>   Trophic state    : Eutrophic
#> 
#> ========================================

Full pipeline in one call

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)
#> ========================================
#>   okBATHTUB Water Quality Summary
#> ========================================
#> 
#>   Segment      : main
#>   Coefficients : walker
#>   Pipeline     : tsi
#> 
#>   -- Hydraulics --
#>   Inflow           : 4.500e+07 m3/yr
#>   Surface area     : 890.0 ha
#>   Mean depth       : 4.20 m
#>   Residence time   : 0.831 yr
#>   Areal water load : 5.06 m/yr
#> 
#>   -- Nutrient Retention --
#>   TP retention     : 0.632  (walker_model1)
#>   TN retention     : 0.553  (walker_model1)
#> 
#>   -- In-Lake Predictions --
#>   TP               : 44.2 ug/L
#>   TN               : 803.8 ug/L
#>   Chlorophyll-a    : 17.68 ug/L
#>   Secchi depth     : 1.06 m
#> 
#>   -- Carlson Trophic State Index --
#>   TSI(TP)          : 58.8
#>   TSI(Chl-a)       : 58.8
#>   TSI(Secchi)      : 59.1
#>   TSI(mean)        : 58.9  (n = 3 components)
#>   Trophic state    : Eutrophic
#> 
#> ========================================

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:

res <- ok_reservoir("Arcadia Lake", exact = TRUE)
res[, c("lake_name", "surface_area_ha", "mean_depth_m", "eco_l3_name",
        "data_quality")]
#>      lake_name surface_area_ha mean_depth_m   eco_l3_name data_quality
#> 1 Arcadia Lake             716          4.6 Cross Timbers            A

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")
#> Trophic state: Eutrophic
cat("Mean TSI:     ", round(result$data$tsi_mean, 1), "\n")
#> Mean TSI:      58.8

For a dataset overview by ecoregion:

ok_reservoir_summary()
#> ========================================
#>   okBATHTUB - ok_reservoirs Coverage
#> ========================================
#> 
#>   Total lakes: 40
#>   Quality A (authoritative source): 39
#>   Quality B (depth estimated):      1
#> 
#>                    eco_l3_name n_lakes n_quality_a n_quality_b area_min_ha
#>                  Cross Timbers      14          14           0         440
#>                Ozark Highlands       7           7           0         440
#>                Arkansas Valley       6           5           1         212
#>  Central Oklahoma/Texas Plains       6           6           0         720
#>             Ouachita Mountains       6           6           0        1620
#>        Southwestern Tablelands       1           1           0        2090
#>  area_max_ha depth_min_m depth_max_m
#>        10570         2.4        11.5
#>        18700         3.8        11.3
#>        41280         3.4         7.0
#>        35840         4.3        11.6
#>         5710         4.4        15.5
#>         2090         3.5         3.5

Multiple tributaries

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

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")
#> FW mean TP inflow: 120 ug/L

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:

obs <- ok_tsi_observed(
  tp_ugl    = 85,
  chla_ugl  = 22,
  secchi_m  = 0.8
)

cat("TSI(TP):    ", round(obs$tsi_tp,     1), "\n")
#> TSI(TP):     68.2
cat("TSI(Chl-a): ", round(obs$tsi_chla,   1), "\n")
#> TSI(Chl-a):  60.9
cat("TSI(Secchi):", round(obs$tsi_secchi, 1), "\n")
#> TSI(Secchi): 63.2
cat("Mean TSI:   ", round(obs$tsi_mean,   1), "\n")
#> Mean TSI:    64.1
cat("Class:      ", obs$trophic_state,         "\n")
#> Class:       Eutrophic

Comparing coefficient sets

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

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)
#>                     scenario tp_inlake  chla tsi_mean
#>               Walker Model 1      44.2 17.68     58.9
#>  Vollenweider/Larsen-Mercier      62.8 29.45     63.4
#>     Oklahoma (Cross Timbers)      44.2 19.83     62.3

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.