Oklahoma Reservoir Workflow

Overview

This vignette demonstrates Oklahoma-focused workflows with okBATHTUB:

Looking up reservoir morphometry from the bundled dataset
Selecting appropriate coefficient sets for different ecoregions
Comparing model predictions to observed monitoring data
Multi-segment reservoir modelling
Load reduction scenario analysis

library(okBATHTUB)

Reservoir lookup

The bundled ok_reservoirs dataset includes morphometric data for major Oklahoma reservoirs compiled from publicly available sources (USACE Tulsa District design memoranda, U.S. Bureau of Reclamation design data, the National Inventory of Dams, and published OWRB reports).

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

Look up any reservoir by name (partial or exact match):

# Partial match
ok_reservoir("Tenkiller")[, c("lake_name", "surface_area_ha", "mean_depth_m",
                              "eco_l3_name")]
#>              lake_name surface_area_ha mean_depth_m        eco_l3_name
#> 4 Tenkiller Ferry Lake            5060         15.5 Ouachita Mountains

# Filter by ecoregion
ct <- ok_reservoir(ecoregion = "Cross Timbers")
nrow(ct)
#> [1] 14
head(ct[, c("lake_name", "surface_area_ha", "mean_depth_m")], 5)
#>           lake_name surface_area_ha mean_depth_m
#> 1      Arcadia Lake             716          4.6
#> 5  Lake Thunderbird            2386          5.5
#> 7     Keystone Lake           10570          7.7
#> 11      Lake Hefner            1015          8.0
#> 12  Lake Overholser             615          2.4

Choosing a coefficient set

Three coefficient families are supported:

Set	Retention model	Chl-a / Secchi
`"walker"`	Walker BATHTUB Model 1	Walker (1985) national
`"vollenweider"`	Vollenweider/Larsen-Mercier	Walker (1985) national
`"oklahoma"`	Walker BATHTUB Model 1	Oklahoma ecoregion-specific

For Oklahoma reservoirs, "oklahoma" is generally preferred because the Chl-a / Secchi regressions are calibrated to Oklahoma data. The retention model is still Walker Model 1, which is the canonical default in the BATHTUB program and is calibrated to U.S. Army Corps of Engineers reservoir data — many Oklahoma reservoirs are Corps projects, so this is appropriate.

If you have reason to prefer a simpler retention model (no ortho-P / total-P data, very rough screening, or comparison to legacy results), "vollenweider" is available. Be aware that this form is not calibrated to Corps of Engineers reservoir data.

Worked example: Arcadia Lake

Arcadia Lake is a 716 ha Cross Timbers reservoir managed by the City of Edmond for water supply. Suppose we have estimated tributary loads from a watershed model:

arcadia <- ok_reservoir("Arcadia Lake", exact = TRUE)

result <- ok_load(
    inflow_m3yr   = 45e6,    # ~36,500 acre-feet/yr
    tp_inflow_ugl = 120,     # flow-weighted mean TP
    tn_inflow_ugl = 1800,    # flow-weighted mean TN
    coefficients  = "oklahoma",
    ecoregion     = "Cross Timbers"
  ) |>
  ok_hydraulics(
    surface_area_ha = arcadia$surface_area_ha,
    mean_depth_m    = arcadia$mean_depth_m
  ) |>
  ok_retention() |>
  ok_inlake()    |>
  ok_tsi()

summary(result)
#> ========================================
#>   okBATHTUB Water Quality Summary
#> ========================================
#> 
#>   Segment      : main
#>   Coefficients : oklahoma
#>   Ecoregion    : Cross Timbers
#>   Pipeline     : tsi
#> 
#>   -- Hydraulics --
#>   Inflow           : 4.500e+07 m3/yr
#>   Surface area     : 716.0 ha
#>   Mean depth       : 4.60 m
#>   Residence time   : 0.732 yr
#>   Areal water load : 6.28 m/yr
#> 
#>   -- Nutrient Retention --
#>   TP retention     : 0.636  (walker_model1)
#>   TN retention     : 0.553  (walker_model1)
#> 
#>   -- In-Lake Predictions --
#>   TP               : 43.7 ug/L
#>   TN               : 805.1 ug/L
#>   Chlorophyll-a    : 19.71 ug/L
#>   Secchi depth     : 0.57 m
#> 
#>   -- Carlson Trophic State Index --
#>   TSI(TP)          : 58.6
#>   TSI(Chl-a)       : 59.8
#>   TSI(Secchi)      : 68.1
#>   TSI(mean)        : 62.2  (n = 3 components)
#>   Trophic state    : Eutrophic
#> 
#> ========================================

Comparing to observed monitoring data

If you have observed in-lake monitoring data, use ok_tsi_observed() to compute observed TSI and compare to predictions:

# Suppose lake monitoring data summary for Arcadia gives these
# growing-season means:
obs <- ok_tsi_observed(
  tp_ugl   = 56,
  chla_ugl = 18,
  secchi_m = 0.95
)

cat("--- Predicted (okBATHTUB) vs Observed ---\n")
#> --- Predicted (okBATHTUB) vs Observed ---
cat(sprintf("TP    : pred = %5.1f, obs = %5.1f ug/L\n",
            result$data$tp_inlake_ugl, 56))
#> TP    : pred =  43.7, obs =  56.0 ug/L
cat(sprintf("Chl-a : pred = %5.2f, obs = %5.2f ug/L\n",
            result$data$chla_ugl, 18))
#> Chl-a : pred = 19.71, obs = 18.00 ug/L
cat(sprintf("Secchi: pred = %5.2f, obs = %5.2f m\n",
            result$data$secchi_m, 0.95))
#> Secchi: pred =  0.57, obs =  0.95 m
cat(sprintf("TSI   : pred = %5.1f, obs = %5.1f\n",
            result$data$tsi_mean, obs$tsi_mean))
#> TSI   : pred =  62.2, obs =  60.6

The predicted-vs-observed comparison is the most important diagnostic for any empirical model. Persistent biases in one direction may indicate that the model coefficients need site-specific calibration: pass a custom named list as coefficients to ok_load() if you have sufficient data to fit your own.

Multi-segment reservoir modelling

Reservoirs with strong longitudinal gradients (Tenkiller, Eufaula, and Texoma are classic Oklahoma examples) are often best modelled as a series of well-mixed segments. Each segment’s outflow feeds the next segment’s inflow:

# A simple 3-segment Tenkiller-like reservoir
result <- ok_segment_chain(
  inflow_m3yr   = 1e9,
  tp_inflow_ugl = 110,
  tn_inflow_ugl = 1700,
  segments = list(
    list(label = "riverine",    surface_area_ha = 1200, mean_depth_m = 4),
    list(label = "transition",  surface_area_ha = 1800, mean_depth_m = 9),
    list(label = "lacustrine",  surface_area_ha = 2060, mean_depth_m = 25)
  ),
  coefficients = "oklahoma",
  ecoregion    = "Ouachita Mountains"
)
#> okBATHTUB: no ecoregion-specific Chl-a coefficients for 'Ouachita Mountains'; using Oklahoma statewide pooled.
#> okBATHTUB: no ecoregion-specific Chl-a coefficients for 'Ouachita Mountains'; using Oklahoma statewide pooled.
#> okBATHTUB: no ecoregion-specific Chl-a coefficients for 'Ouachita Mountains'; using Oklahoma statewide pooled.

ok_segment_summary(result)
#>      segment tp_inflow_ugl tp_inlake_ugl tp_retention tn_inlake_ugl  chla_ugl
#> 1   riverine     110.00000      72.75332        0.339     1341.9538 25.161363
#> 2 transition      72.75332      38.98392        0.464      862.1997 16.549299
#> 3 lacustrine      38.98392      17.62827        0.548      451.2492  9.712472
#>    secchi_m   tsi_tp tsi_chla tsi_secchi tsi_mean trophic_state hrt_yr
#> 1 0.5326092 65.96961 62.24029   69.07783 65.76258     Eutrophic  0.048
#> 2 0.6658716 56.97261 58.13023   65.85995 60.32093     Eutrophic  0.162
#> 3 0.8846128 45.52824 52.90216   61.76674 53.39905     Eutrophic  0.515
#>   surface_area_ha mean_depth_m
#> 1            1200            4
#> 2            1800            9
#> 3            2060           25

In-lake TP decreases moving down-reservoir as nutrients settle, which is the expected pattern for an unstratified-to-stratified reservoir gradient.

Load reduction scenarios

Use ok_scenario_sweep() to evaluate a range of nutrient load reductions and identify the reduction needed to meet a TSI target:

baseline <- ok_load(
    inflow_m3yr   = 45e6,
    tp_inflow_ugl = 120,
    coefficients  = "oklahoma",
    ecoregion     = "Cross Timbers"
  ) |>
  ok_hydraulics(surface_area_ha = 716, mean_depth_m = 4.6)

sweep <- ok_scenario_sweep(
  baseline,
  max_reduction_pct = 60,
  step_pct          = 10
)

# Print key columns
print(sweep[, c("tp_reduction_pct", "tp_inflow_ugl",
                "tp_inlake_ugl", "chla_ugl", "tsi_mean",
                "trophic_state")],
      row.names = FALSE)
#>  tp_reduction_pct tp_inflow_ugl tp_inlake_ugl chla_ugl tsi_mean trophic_state
#>                 0           120          43.7    19.71     62.2     Eutrophic
#>                10           108          41.0    18.94     61.7     Eutrophic
#>                20            96          38.1    18.10     61.0     Eutrophic
#>                30            84          35.0    17.19     60.3     Eutrophic
#>                40            72          31.7    16.18     59.5     Eutrophic
#>                50            60          28.2    15.04     58.5     Eutrophic
#>                60            48          24.3    13.73     57.3     Eutrophic

References

Carlson, R.E. (1977). A trophic state index for lakes. Limnology and Oceanography, 22(2), 361-369.
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.