OK-HAWQS/SWAT and okBATHTUB are complementary models designed to work together in a two-model nutrient management workflow:
Neither model alone spans the full management question. Together they answer: “If we implement BMPs that reduce watershed TP export by X%, what will the reservoir’s chlorophyll-a and trophic state be?”
Watershed Reservoir
─────────────────────────────────── ──────────────────────────────────
Land use / management inputs Tributary loads from HAWQS
↓ ↓
OK-HAWQS/SWAT ok_load() ──→ ok_hydraulics()
↓ ↓
Daily streamflow + loads ok_retention() ──→ ok_inlake()
↓ ↓
Annual/seasonal load summaries ok_tsi() ──→ ok_plot_response()
↓ ↓
[Handoff point] In-lake TP, Chl-a, Secchi, TSIThe handoff occurs at the watershed outlet / reservoir inlet. OK-HAWQS produces annual or seasonal mean tributary concentrations and flow volumes; okBATHTUB takes these as inputs.
OK-HAWQS/SWAT produces output files in the TxtInOut
directory. The relevant files for reservoir modelling are:
| File | Contents | okBATHTUB use |
|---|---|---|
output.rch |
Streamflow and constituent loads at each reach | Inflow volume and TP/TN loads |
output.sub |
Subbasin water balance | Watershed area confirmation |
output.hru |
HRU-level outputs | BMP scenario comparison |
The output.rch file uses fixed-width format with columns
including FLOW_OUTcms (flow in m³/s), NO3_OUT
(nitrate), ORGN_OUT (organic N), SOLP_OUT
(soluble P), and ORGP_OUT (organic P).
# Read SWAT output.rch file
# (Adjust column positions for your SWAT version)
read_swat_rch <- function(rch_path) {
# Read fixed-width format (requires the 'readr' package)
raw <- readr::read_fwf(
rch_path,
readr::fwf_cols(
rch_id = c(1, 4),
year = c(5, 9),
month = c(10, 12),
day = c(13, 16),
area_km2 = c(17, 26),
flow_cms = c(27, 38), # FLOW_OUTcms
no3_kgha = c(83, 94), # NO3_OUT
orgn_kgha = c(95, 106), # ORGN_OUT
solp_kgha = c(107, 118), # SOLP_OUT
orgp_kgha = c(119, 130) # ORGP_OUT
),
skip = 9,
col_types = "iiiiidddddd"
)
raw
}This example demonstrates the analytical approach typical of state-level Clean Water Act Section 604(b) watershed modeling projects. We use simulated OK-HAWQS output representing annual TP and flow loads at the reservoir inlet for a hypothetical sensitive water supply reservoir.
# Simulated OK-HAWQS annual output at reservoir inlet
# Replace with actual output.rch data in production use
hawqs_annual <- data.frame(
year = 2010:2022,
flow_m3yr = c(38.2, 52.1, 29.4, 44.7, 61.3, 35.8, 41.9,
48.2, 33.6, 57.4, 42.1, 36.8, 49.3) * 1e6,
tp_load_kgyr = c(4584, 7314, 2823, 5811, 9195, 3938, 5447,
6253, 3427, 8296, 5473, 3938, 6657),
tn_load_kgyr = c(69120, 98990, 47100, 84780, 127280, 57890, 79860,
91580, 53420, 114800, 80020, 58920, 95280)
)
# Compute flow-weighted mean concentrations
# Concentration (µg/L) = load (kg/yr) / volume (m³/yr) * 1e6
hawqs_annual$tp_conc_ugl <- hawqs_annual$tp_load_kgyr / hawqs_annual$flow_m3yr * 1e6
hawqs_annual$tn_conc_ugl <- hawqs_annual$tn_load_kgyr / hawqs_annual$flow_m3yr * 1e6
cat("--- OK-HAWQS Annual Load Summary ---\n")
#> --- OK-HAWQS Annual Load Summary ---
cat(sprintf("Mean annual flow: %.1f million m³/yr\n",
mean(hawqs_annual$flow_m3yr) / 1e6))
#> Mean annual flow: 43.9 million m³/yr
cat(sprintf("Mean annual TP load: %.0f kg/yr\n",
mean(hawqs_annual$tp_load_kgyr)))
#> Mean annual TP load: 5627 kg/yr
cat(sprintf("FWM TP conc: %.1f µg/L\n",
mean(hawqs_annual$tp_conc_ugl)))
#> FWM TP conc: 125.0 µg/L
cat(sprintf("FWM TN conc: %.1f µg/L\n",
mean(hawqs_annual$tn_conc_ugl)))
#> FWM TN conc: 1825.5 µg/LThe key conversion is from HAWQS load (kg/yr) and flow (m³/yr) to okBATHTUB inflow concentration (µg/L):
\[C_{inflow} (\mu g/L) = \frac{L (kg/yr) \times 10^6}{Q (m^3/yr)}\]
# Use long-term mean conditions for steady-state BATHTUB
mean_flow_m3yr <- mean(hawqs_annual$flow_m3yr)
mean_tp_ugl <- mean(hawqs_annual$tp_conc_ugl)
mean_tn_ugl <- mean(hawqs_annual$tn_conc_ugl)
cat(sprintf("BATHTUB inputs:\n"))
#> BATHTUB inputs:
cat(sprintf(" Inflow volume: %.2f million m³/yr\n", mean_flow_m3yr / 1e6))
#> Inflow volume: 43.91 million m³/yr
cat(sprintf(" TP inflow: %.1f µg/L\n", mean_tp_ugl))
#> TP inflow: 125.0 µg/L
cat(sprintf(" TN inflow: %.1f µg/L\n", mean_tn_ugl))
#> TN inflow: 1825.5 µg/L# Reservoir morphometry: a representative Cross Timbers reservoir
# surface_area_ha = 890, mean_depth_m = 4.2 (hardcoded for reproducibility)
baseline <- ok_load(
inflow_m3yr = mean_flow_m3yr,
tp_inflow_ugl = mean_tp_ugl,
tn_inflow_ugl = mean_tn_ugl,
coefficients = "oklahoma",
ecoregion = "Cross Timbers",
segment_label = "baseline"
) |>
ok_hydraulics(
surface_area_ha = 890,
mean_depth_m = 4.2
)
result_baseline <- baseline |>
ok_retention() |>
ok_inlake() |>
ok_tsi()
summary(result_baseline)
#> ========================================
#> okBATHTUB Water Quality Summary
#> ========================================
#>
#> Segment : baseline
#> Coefficients : oklahoma
#> Ecoregion : Cross Timbers
#> Pipeline : tsi
#>
#> -- Hydraulics --
#> Inflow : 4.391e+07 m3/yr
#> Surface area : 890.0 ha
#> Mean depth : 4.20 m
#> Residence time : 0.851 yr
#> Areal water load : 4.93 m/yr
#>
#> -- Nutrient Retention --
#> TP retention : 0.639 (walker_model1)
#> TN retention : 0.557 (walker_model1)
#>
#> -- In-Lake Predictions --
#> TP : 45.2 ug/L
#> TN : 808.2 ug/L
#> Chlorophyll-a : 20.11 ug/L
#> Secchi depth : 0.56 m
#>
#> -- Carlson Trophic State Index --
#> TSI(TP) : 59.1
#> TSI(Chl-a) : 60.0
#> TSI(Secchi) : 68.3
#> TSI(mean) : 62.5 (n = 3 components)
#> Trophic state : Eutrophic
#>
#> ========================================OK-HAWQS/SWAT is designed to simulate the effect of conservation practices (cover crops, reduced tillage, buffer strips, wetland restoration) on watershed loading. Each HAWQS scenario produces a modified load — which feeds directly into okBATHTUB.
Here we demonstrate using simulated HAWQS BMP scenario outputs:
# Simulated HAWQS BMP scenario outputs
# Each represents a different watershed management alternative
hawqs_scenarios <- list(
list(
label = "Baseline (current conditions)",
flow_m3yr = mean_flow_m3yr,
tp_load_kgyr = mean(hawqs_annual$tp_load_kgyr),
tn_load_kgyr = mean(hawqs_annual$tn_load_kgyr)
),
list(
label = "10% cropland cover crops",
flow_m3yr = mean_flow_m3yr * 0.97, # slight flow reduction
tp_load_kgyr = mean(hawqs_annual$tp_load_kgyr) * 0.88,
tn_load_kgyr = mean(hawqs_annual$tn_load_kgyr) * 0.92
),
list(
label = "30% cropland cover crops + buffer strips",
flow_m3yr = mean_flow_m3yr * 0.94,
tp_load_kgyr = mean(hawqs_annual$tp_load_kgyr) * 0.72,
tn_load_kgyr = mean(hawqs_annual$tn_load_kgyr) * 0.78
),
list(
label = "Full BMP suite (TMDL alternative)",
flow_m3yr = mean_flow_m3yr * 0.91,
tp_load_kgyr = mean(hawqs_annual$tp_load_kgyr) * 0.55,
tn_load_kgyr = mean(hawqs_annual$tn_load_kgyr) * 0.62
)
)
# Convert each HAWQS scenario to okBATHTUB concentrations and run pipeline
scenario_results <- lapply(hawqs_scenarios, function(sc) {
tp_ugl <- sc$tp_load_kgyr / sc$flow_m3yr * 1e6
tn_ugl <- sc$tn_load_kgyr / sc$flow_m3yr * 1e6
r <- ok_load(
inflow_m3yr = sc$flow_m3yr,
tp_inflow_ugl = tp_ugl,
tn_inflow_ugl = tn_ugl,
coefficients = "oklahoma",
ecoregion = "Cross Timbers"
) |>
ok_hydraulics(
surface_area_ha = 890,
mean_depth_m = 4.2
) |>
ok_retention() |>
ok_inlake() |>
ok_tsi()
data.frame(
scenario = sc$label,
tp_inflow_ugl = round(tp_ugl, 1),
tp_reduction_pct = round(100 * (1 - tp_ugl / mean_tp_ugl), 1),
tp_inlake_ugl = round(r$data$tp_inlake_ugl, 1),
chla_ugl = round(r$data$chla_ugl, 2),
secchi_m = round(r$data$secchi_m, 2),
tsi_mean = round(r$data$tsi_mean, 1),
trophic_state = r$data$trophic_state,
stringsAsFactors = FALSE
)
})
scenario_df <- do.call(rbind, scenario_results)
print(scenario_df)
#> scenario tp_inflow_ugl tp_reduction_pct
#> 1 Baseline (current conditions) 128.2 -2.6
#> 2 10% cropland cover crops 116.3 7.0
#> 3 30% cropland cover crops + buffer strips 98.2 21.5
#> 4 Full BMP suite (TMDL alternative) 77.5 38.0
#> tp_inlake_ugl chla_ugl secchi_m tsi_mean trophic_state
#> 1 45.8 20.30 0.56 62.6 Eutrophic
#> 2 43.1 19.53 0.57 62.1 Eutrophic
#> 3 38.7 18.27 0.59 61.2 Eutrophic
#> 4 33.2 16.64 0.62 59.9 EutrophicOklahoma’s flood-dominated hydrology means that annual TP loading varies substantially between wet and dry years. A robust BATHTUB analysis should bracket this variability rather than relying on a single mean year:
# Run BATHTUB for each year in the HAWQS record
annual_list <- lapply(seq_len(nrow(hawqs_annual)), function(i) {
yr <- hawqs_annual[i, ]
tp_ugl <- yr$tp_conc_ugl
tn_ugl <- yr$tn_conc_ugl
r <- ok_load(
inflow_m3yr = yr$flow_m3yr,
tp_inflow_ugl = tp_ugl,
tn_inflow_ugl = tn_ugl,
coefficients = "oklahoma",
ecoregion = "Cross Timbers"
) |>
ok_hydraulics(
surface_area_ha = 890,
mean_depth_m = 4.2
) |>
ok_retention() |>
ok_inlake() |>
ok_tsi()
data.frame(
year = yr$year,
flow_m3yr = yr$flow_m3yr / 1e6,
tp_inflow_ugl = round(tp_ugl, 1),
tp_inlake_ugl = round(r$data$tp_inlake_ugl, 1),
chla_ugl = round(r$data$chla_ugl, 2),
secchi_m = round(r$data$secchi_m, 2),
tsi_mean = round(r$data$tsi_mean, 1),
trophic_state = r$data$trophic_state
)
})
annual_results <- do.call(rbind, annual_list)
cat("--- Interannual Range ---\n")
#> --- Interannual Range ---
cat(sprintf("TSI range: %.1f - %.1f\n",
min(annual_results$tsi_mean), max(annual_results$tsi_mean)))
#> TSI range: 60.3 - 63.7
cat(sprintf("Chl-a range: %.2f - %.2f µg/L\n",
min(annual_results$chla_ugl), max(annual_results$chla_ugl)))
#> Chl-a range: 17.17 - 22.07 µg/L
cat(sprintf("Secchi range: %.2f - %.2f m\n",
min(annual_results$secchi_m), max(annual_results$secchi_m)))
#> Secchi range: 0.54 - 0.61 mThe load-response curve places the HAWQS scenarios in the context of the full TP range:
ok_plot_response(
baseline,
response = "tsi",
target_class = "mesotrophic",
current_tp = mean_tp_ugl,
lake_name = "Cross Timbers Reservoir -- OK-HAWQS/okBATHTUB Linkage"
)
Temporal averaging — BATHTUB is a steady-state model. Use the growing-season (May–October) mean for response variables and the water-year mean for loads. For interannual analysis, run BATHTUB separately for each year in the HAWQS record.
Wet vs dry year stratification — consider running separate BATHTUB analyses for wet-year (above median flow) and dry-year (below median flow) conditions from your HAWQS record to bracket uncertainty.
Load vs concentration — HAWQS produces loads
(kg/yr). Convert to concentration (µg/L) using
C = L / Q * 1e6 before calling ok_load().
Verify that the resulting concentration is physically plausible for the
watershed — very high concentrations with low flows may indicate
numerical artefacts in the HAWQS calibration.
Retention sensitivity — Walker BATHTUB Model 1 retention is sensitive to hydraulic residence time, which depends on inflow volume. In drought years, lower inflow → longer residence time → higher retention → lower in-lake TP. Check that this behaviour is captured when comparing wet and dry year scenarios.
Uncertainty propagation — both HAWQS and BATHTUB have calibration uncertainty. Report predictions as ranges (e.g., ± one standard deviation across years) rather than single values, particularly for regulatory applications.