--- title: "Linking OK-HAWQS/SWAT Outputs to okBATHTUB" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Linking OK-HAWQS/SWAT Outputs to okBATHTUB} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, dpi = 120 ) ``` ## Overview OK-HAWQS/SWAT and okBATHTUB are complementary models designed to work together in a two-model nutrient management workflow: - **OK-HAWQS/SWAT** simulates watershed hydrology and nutrient loading — translating land management decisions into tributary flow, sediment, and nutrient yields at the watershed outlet - **okBATHTUB** simulates in-lake response — translating those tributary loads into predicted reservoir water quality 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?"* ```{r load} library(okBATHTUB) ``` ## The Two-Model Workflow ``` 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, TSI ``` The 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 Output File Structure 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). ## Reading OK-HAWQS Output ```{r read_hawqs, eval=FALSE} # 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 } ``` ## Simulated Example: 604(b) Sensitive Water Supply Workflow 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. ```{r simulate_hawqs} # 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") cat(sprintf("Mean annual flow: %.1f million m³/yr\n", mean(hawqs_annual$flow_m3yr) / 1e6)) cat(sprintf("Mean annual TP load: %.0f kg/yr\n", mean(hawqs_annual$tp_load_kgyr))) cat(sprintf("FWM TP conc: %.1f µg/L\n", mean(hawqs_annual$tp_conc_ugl))) cat(sprintf("FWM TN conc: %.1f µg/L\n", mean(hawqs_annual$tn_conc_ugl))) ``` ## Converting HAWQS Loads to okBATHTUB Inputs The 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)}$$ ```{r hawqs_to_bathtub} # 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")) cat(sprintf(" Inflow volume: %.2f million m³/yr\n", mean_flow_m3yr / 1e6)) cat(sprintf(" TP inflow: %.1f µg/L\n", mean_tp_ugl)) cat(sprintf(" TN inflow: %.1f µg/L\n", mean_tn_ugl)) ``` ## Baseline BATHTUB Run ```{r baseline_run} # 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) ``` ## BMP Scenario Analysis 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: ```{r bmp_scenarios} # 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) ``` ## Interannual Variability Analysis Oklahoma'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: ```{r interannual} # 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") cat(sprintf("TSI range: %.1f - %.1f\n", min(annual_results$tsi_mean), max(annual_results$tsi_mean))) cat(sprintf("Chl-a range: %.2f - %.2f µg/L\n", min(annual_results$chla_ugl), max(annual_results$chla_ugl))) cat(sprintf("Secchi range: %.2f - %.2f m\n", min(annual_results$secchi_m), max(annual_results$secchi_m))) ``` ## Load-Response Curve with HAWQS Context The load-response curve places the HAWQS scenarios in the context of the full TP range: ```{r response_plot, fig.height=5.5, eval=requireNamespace("ggplot2", quietly=TRUE) && exists("baseline")} ok_plot_response( baseline, response = "tsi", target_class = "mesotrophic", current_tp = mean_tp_ugl, lake_name = "Cross Timbers Reservoir -- OK-HAWQS/okBATHTUB Linkage" ) ``` ## Best Practices for HAWQS-BATHTUB 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. ## References - Soil and Water Assessment Tool (SWAT) documentation: https://swat.tamu.edu - OK-HAWQS: Oklahoma application of HAWQS 2.0 (USGS/EPA national platform) - Walker, W.W. (1996). *Simplified Procedures for Eutrophication Assessment and Prediction: User Manual*. U.S. Army Corps of Engineers, Instruction Report W-96-2. - U.S. EPA Clean Water Act Section 604(b) Water Quality Management Planning Grant program: https://www.epa.gov/nps/water-quality-management-planning-grants