Getting Started with okBATHTUB

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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.

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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.