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.

Linking OK-HAWQS/SWAT Outputs to okBATHTUB

Overview

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?”

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

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

# 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/L

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)}\]

# 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

Baseline BATHTUB 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)
#> ========================================
#>   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
#> 
#> ========================================

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:

# 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     Eutrophic

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:

# 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 m

Load-Response Curve with HAWQS Context

The 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"
)

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

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.