Calibration Frameworks

OpenWQ provides a comprehensive Python calibration framework for parameter optimization and sensitivity analysis. The framework is designed for computationally expensive models, supporting both local (Docker) and HPC (Apptainer/SLURM) execution environments.

Key features:

DDS Optimization: Dynamically Dimensioned Search algorithm, efficient for expensive models (200-1000 evaluations)
Sensitivity Analysis: Morris screening and Sobol variance-based methods for parameter importance ranking
All Parameter Types: Supports BGC rates, PHREEQC geochemistry, sorption isotherms, sediment transport, source/sink loads, and more
HPC Support: SLURM/PBS job arrays for parallel evaluation on clusters
Checkpointing: Save and resume long-running calibrations

Location

The calibration framework is located at:

openwq/supporting_scripts/Calibration/

Key files:

calibration_config_template.py - User configuration template (copy and edit)
calibration_lib/calibration_driver.py - Main entry point
calibration_lib/ - Core library modules (including observation_data/)

Quick Start

Copy the configuration template:

cd openwq/supporting_scripts/Calibration/
cp calibration_config_template.py my_calibration.py

Edit paths and parameters in my_calibration.py
Run calibration:
```
python my_calibration.py
```
Resume if interrupted:
```
python my_calibration.py --resume
```

Run sensitivity analysis only:

python my_calibration.py --sensitivity-only

Validate configuration (dry run):
```
python my_calibration.py --dry-run
```

Configuration Overview

The configuration file has seven main sections:

Path Configuration - Model config, observations, working directory
Container Runtime - Docker (local) or Apptainer (HPC) settings
Calibration Parameters - Parameters to optimize with bounds
Calibration Settings - Algorithm, max evaluations, objective function
Sensitivity Analysis - Morris/Sobol settings
HPC Settings - SLURM/PBS job configuration
Run Calibration - Entry point

Supported Parameter Types

The framework supports calibration of all major OpenWQ parameter types:

BGC Parameters (NATIVE_BGC_FLEX)

{
    "name": "k_nitrification",
    "file_type": "bgc_json",
    "path": ["CYCLING_FRAMEWORKS", "N_cycle", "3", "parameter_values", "k"],
    "initial": 0.03,
    "bounds": (0.001, 0.5),
    "transform": "log"
}

PHREEQC Geochemistry Parameters

# Initial solution concentrations
{
    "name": "initial_NO3",
    "file_type": "phreeqc_pqi",
    "path": {"block": "SOLUTION", "species": "N(5)"},
    "initial": 2.0,
    "bounds": (0.05, 50.0)
}

# Equilibrium phase partial pressures
{
    "name": "pCO2_log",
    "file_type": "phreeqc_pqi",
    "path": {"block": "EQUILIBRIUM_PHASES", "phase": "CO2(g)", "field": "si"},
    "initial": -3.5,
    "bounds": (-4.5, -2.5)
}

Sorption Isotherm Parameters

# Freundlich isotherm
{
    "name": "NH4_Kfr",
    "file_type": "sorption_json",
    "path": {"module": "FREUNDLICH", "species": "NH4-N", "param": "Kfr"},
    "initial": 0.5,
    "bounds": (0.01, 10.0),
    "transform": "log"
}

# Langmuir isotherm
{
    "name": "NH4_qmax",
    "file_type": "sorption_json",
    "path": {"module": "LANGMUIR", "species": "NH4-N", "param": "qmax_mg_per_kg"},
    "initial": 100.0,
    "bounds": (10.0, 500.0)
}

Sediment Transport Parameters

# HYPE_HBVSED module
{
    "name": "erosion_index",
    "file_type": "sediment_json",
    "path": {"module": "HYPE_HBVSED", "param": "EROSION_INDEX"},
    "initial": 0.4,
    "bounds": (0.05, 5.0)
}

# HYPE_MMF module
{
    "name": "cohesion",
    "file_type": "sediment_json",
    "path": {"module": "HYPE_MMF", "param": "COHESION"},
    "initial": 7.5,
    "bounds": (1.0, 20.0)
}

Source/Sink Parameters

# CSV load scaling
{
    "name": "fertilizer_N_scale",
    "file_type": "ss_csv_scale",
    "path": {"species": "all"},
    "initial": 1.0,
    "bounds": (0.1, 5.0)
}

# Copernicus LULC export coefficients
{
    "name": "cropland_TN_coeff",
    "file_type": "ss_copernicus_static",
    "path": {"lulc_class": 10, "species": "TN"},
    "initial": 15.0,
    "bounds": (1.0, 30.0)
}

# Climate response parameters
{
    "name": "precip_scaling_power",
    "file_type": "ss_climate_param",
    "path": "ss_climate_precip_scaling_power",
    "initial": 1.0,
    "bounds": (0.5, 2.0)
}

Transport Parameters

{
    "name": "dispersion_x",
    "file_type": "transport_json",
    "path": {"param": "DISPERSION_X_M2_S"},
    "initial": 0.3,
    "bounds": (0.01, 10.0),
    "transform": "log"
}

Parameter Options

Each parameter definition supports these options:

Option	Description
`name`	Unique identifier (used in results)
`file_type`	Type of configuration file (see list above)
`path`	Location of parameter within the file
`initial`	Starting value for optimization
`bounds`	(min, max) tuple defining search space
`transform`	`"linear"` (default) or `"log"` for log-space optimization
`dtype`	`"float"` (default) or `"int"` for integer parameters

Objective Functions

The framework supports four objective functions:

RMSE - Root Mean Square Error
NSE - Nash-Sutcliffe Efficiency (1 - NSE for minimization)
KGE - Kling-Gupta Efficiency (1 - KGE for minimization) - Recommended
PBIAS - Percent Bias (absolute value for minimization)

Configure in your calibration file:

objective_function = "KGE"

# Multi-species weighting
objective_weights = {
    "NO3-N": 1.0,
    "NH4-N": 0.5,
}

# Target species and locations
calibration_targets = {
    "species": ["NO3-N", "NH4-N"],
    "reach_ids": "all",  # or [1200014181, 200014181]
    "compartments": ["RIVER_NETWORK_REACHES"]
}

Temporal Resolution

The calibration framework supports flexible temporal resolution for objective function calculation and performance evaluation. Both model outputs and observations are aggregated to the specified resolution before computing metrics.

Available Resolutions:

native - Use original model/observation timestamps (no aggregation)
daily - Aggregate to daily averages
weekly - Aggregate to weekly averages
monthly - Aggregate to monthly averages
yearly - Aggregate to yearly averages

Aggregation Methods:

mean - Average value (default, recommended for concentrations)
sum - Total sum (useful for loads/fluxes)
median - Median value (robust to outliers)
min / max - Minimum/maximum values

Configuration:

# Temporal resolution for calibration
temporal_resolution = "monthly"    # Options: native, daily, weekly, monthly, yearly
aggregation_method = "mean"        # Options: mean, sum, median, min, max

When to Use Different Resolutions:

Resolution	Recommended Use Case
`native`	High-frequency data, sub-daily dynamics important
`daily`	Standard choice when daily patterns matter
`weekly`	Smooth out day-to-day variability
`monthly`	Seasonal patterns, reduce noise from sparse observations
`yearly`	Long-term trends, annual budgets

How Aggregation Works:

Observations are grouped by reach_id, species, and temporal period
Model outputs are extracted and matched to observation timestamps
Both are aggregated to the specified resolution using the chosen method
Objective functions are computed on the aggregated data
Performance plots show both native and aggregated comparisons

Example Use Case:

If you have sparse monthly grab samples but an hourly model, using temporal_resolution = "monthly" ensures fair comparison:

# Monthly calibration for sparse observations
temporal_resolution = "monthly"
aggregation_method = "mean"

# This aggregates both model outputs and observations to monthly means
# before computing KGE, NSE, RMSE, etc.

Performance Plots:

When temporal aggregation is enabled, the calibration framework generates additional diagnostic plots:

Time series comparison (observed vs simulated at specified resolution)
Scatter plot with 1:1 line and performance metrics
Residual analysis (temporal patterns in errors)

These are saved in results/performance_plots/ after calibration completes.

Observation Data Format

Observation data must be in CSV format with these columns:

datetime,reach_id,species,value,units,source,uncertainty,quality_flag
2010-06-01 00:00:00,1200014181,NO3-N,2.50,mg/l,USGS_station_A,0.25,GOOD
2010-06-01 00:00:00,1200014181,NH4-N,0.15,mg/l,USGS_station_A,0.02,GOOD

Required columns:

datetime - Observation timestamp (YYYY-MM-DD HH:MM:SS)
reach_id - River network reach ID matching model output
species - Chemical species name (must match model species)
value - Measured concentration value
units - Concentration units (e.g., mg/l)

Optional columns:

source - Data source identifier
uncertainty - Measurement uncertainty
quality_flag - Quality flag (GOOD, SUSPECT, BAD)

An example file is provided at observation_data/example_observations.csv.

GRQA Database Integration

The calibration framework can automatically extract observation data from the Global River Water Quality Archive (GRQA) database (https://zenodo.org/records/15335450). This is recommended for large-scale calibrations where manually preparing observation data would be impractical.

Enabling GRQA Extraction

In your calibration configuration file:

# Enable GRQA data extraction
use_grqa = True

# Configure GRQA extraction
grqa_config = {
    # River network shapefile for spatial matching
    "river_network_shapefile": "/path/to/river_network.shp",

    # Column in shapefile containing reach IDs
    "reach_id_column": "seg_id",

    # Bounding box to filter stations [min_lon, min_lat, max_lon, max_lat]
    "bounding_box": [-125, 24, -66, 50],  # e.g., CONUS

    # Time period for observations
    "start_date": "2000-01-01",
    "end_date": "2020-12-31",

    # Maximum distance (meters) to match stations to river reaches
    "max_station_distance_m": 500,

    # Species mapping: GRQA parameter names -> Model species names
    "species_mapping": {
        "NO3": "NO3-N",      # Nitrate
        "NH4": "NH4-N",      # Ammonium
        "TN": "TN",          # Total Nitrogen
        "PO4": "PO4-P",      # Orthophosphate
        "TP": "TP",          # Total Phosphorus
        "DO": "DO",          # Dissolved Oxygen
    },
}

How GRQA Extraction Works

Download: Only required GRQA parameter files are downloaded based on your species mapping
Filter: Stations are filtered by bounding box and time period
Match: Stations are spatially matched to model river reaches using buffer distance
Convert: Units are converted and data is formatted for calibration
Output: Calibration-ready CSV file is generated automatically

Command-Line Options

# Extract GRQA data only (without running calibration)
python my_calibration.py --extract-grqa-only

# Normal run with automatic GRQA extraction
python my_calibration.py

Species Mapping

The species mapping dictionary maps GRQA parameter names to your model species names. Only parameters listed in this mapping will be downloaded and processed. Available GRQA parameters include:

Nitrogen: NO3, NH4, TN, NO2, DON
Phosphorus: PO4, TP, DP
Carbon: DOC, TOC, BOD, COD
Other: DO, TSS, Chl-a, EC, pH, TDS, Temp

Dependencies

GRQA extraction requires additional Python packages:

pip install geopandas requests shapely

Output Files

GRQA extraction creates the following files in your calibration workspace:

grqa_stations.csv - Extracted station metadata
grqa_observations.csv - Raw observations from GRQA
calibration_observations.csv - Formatted for calibration (auto-used)
station_reach_matching.csv - Station-to-reach spatial matches

DDS Optimization Algorithm

The Dynamically Dimensioned Search (DDS) algorithm is a global optimization method designed for computationally expensive models. Key features:

Efficient convergence: Typically requires 200-500 evaluations
Adaptive search: Starts with global exploration, progressively focuses on local refinement
No gradient needed: Uses perturbation-based search

Configure DDS settings:

algorithm = "DDS"
max_evaluations = 500
random_seed = 42  # For reproducibility

Sensitivity Analysis

For calibrations with many parameters (10+), run sensitivity analysis first to identify influential parameters.

Morris Screening (fast, qualitative)

run_sensitivity_first = True
sensitivity_method = "morris"
sensitivity_morris_trajectories = 10
sensitivity_morris_levels = 4
sensitivity_threshold = 0.1  # Relative to max mu_star

Sobol Analysis (slower, quantitative)

sensitivity_method = "sobol"
sensitivity_sobol_samples = 1024

Recommended HPC Workflow

For large-scale calibrations with dozens of parameters:

Morris Screening (~200 evaluations)
- Identifies non-influential parameters
- Typically eliminates 50-70% of parameters
Sobol Sensitivity (~1000 evaluations on reduced set)
- Quantifies parameter importance and interactions
- Further reduces to 10-15 most influential parameters
DDS Optimization (~300-500 evaluations on final set)
- Calibrates only influential parameters
- Fixed non-influential parameters at defaults

Total: ~1500-2000 evaluations (vs 10,000+ without screening)

HPC Configuration (SLURM)

For running on HPC clusters with Apptainer:

# Container runtime
container_runtime = "apptainer"
apptainer_sif_path = "/path/to/openwq.sif"
apptainer_bind_path = "/scratch/user/openwq_code:/code"

# HPC settings
hpc_enabled = True
hpc_scheduler = "slurm"
hpc_partition = "standard"
hpc_walltime = "24:00:00"
hpc_nodes = 1
hpc_tasks_per_node = 4
hpc_memory = "8G"
hpc_max_concurrent_jobs = 50

SLURM job templates are provided in hpc_templates/:

slurm_array_template.sh - For sensitivity analysis (parallel array)
slurm_worker_template.sh - For DDS optimization (worker pool)
pbs_array_template.sh - For PBS/Torque schedulers

Checkpointing and Resume

The framework automatically saves checkpoints after each evaluation. To resume an interrupted calibration:

python my_calibration.py --resume

Checkpoint files are saved in your working directory:

calibration_state.json - Optimization state
calibration_history.pkl - Full evaluation history
rng_state.pkl - Random number generator state

Results and Outputs

After calibration completes, results are saved to results/:

best_parameters.json - Optimal parameter values
calibration_history.json - All evaluated parameter sets and objectives
parameter_definitions.json - Parameter metadata, bounds, and transforms
matched_data.csv - Observation-model matched data pairs per species
calibration_report.html - Interactive HTML report (per-variant)
basin_report.html - Per-basin multi-variant report (auto-generated)
sensitivity_results.json - Sensitivity analysis results (if run)

Interactive HTML Calibration Reports

The calibration framework automatically generates self-contained HTML reports with interactive diagnostics powered by Plotly.js. These reports are produced at the end of each calibration run — no additional steps are needed.

Per-Variant Report (calibration_report.html):

Each completed calibration variant produces a report containing:

6 Interactive Plotly.js Charts: Convergence curve, parameter evolution, time series comparison (observed vs simulated), scatter plot with 1:1 line, residual analysis, and parameter sensitivity ranking
Best Parameters Table: Optimal values alongside initial values, lower/upper bounds, transform type, parameter type, and a visual position-in-range bar
Per-Species Metrics: KGE, NSE, RMSE, PBIAS computed for each target species
Configuration Summary: Algorithm settings, temporal resolution, target species
Dark/Light Theme Toggle: Switch between visual themes
Sidebar Navigation: Quick-jump to any section with scroll-spy highlighting

The report is entirely self-contained (embedded CSS, JavaScript, and data) and can be shared via email or viewed in any modern browser without a web server.

Per-Basin Report (basin_report.html):

When calibrating multiple variants for the same basin (e.g., Variant A: NATIVE_BGC_FLEX, Variant B: SWAT, Variant C: thermodynamic, Variant D: PHREEQC), a consolidated basin report is auto-generated that compares all variants side by side.

Contents:

Overview KPIs: Best objective value per variant, number of evaluations, parameter counts
Variant Comparison Table: Per-species KGE badges color-coded by performance
Per-Variant Detail Cards: Full parameter tables showing Best Value, Initial, Lower Bound, Upper Bound, Transform, Type, and a visual position bar
Interactive Basin Map: Shared spatial context (see below)
Links to Detailed Reports: Click through to each variant’s full report

Auto-Discovery Mechanism: The basin report is generated by scanning sibling workspace directories. Workspaces must follow the naming convention workspace_{basin_id}_{variant} (e.g., workspace_CAN_01AM001_meso_A, workspace_CAN_01AM001_meso_B). The framework detects all completed variants and aggregates their results.

Interactive Basin Maps

Calibration reports include interactive Leaflet.js maps that display the basin’s spatial context. Maps are embedded directly in the HTML reports and require no external server or GIS software.

Map Layers:

HRU Polygons: Hydrologic Response Units colored by area quintiles (indigo palette). Click any polygon to see its attributes (area, HRU ID, elevation)
River Network: Stream segments styled by Strahler order — thicker, darker lines for higher-order streams
Observation Stations: Red circle markers at monitoring locations with popup info (station name, reach ID, available species)

Map Controls:

3 selectable basemaps: CARTO Light, OpenTopoMap, Esri Satellite
Layer toggle (show/hide HRUs, river network, stations)
Legend, scale bar, and zoom controls
Theme-aware: adapts to dark/light mode

Basin Info Grid: Below the map, a summary grid shows: total HRUs, river reaches, total area (km²), network length (km), maximum Strahler order, and number of observation stations.

Data Requirements: Maps are generated from GeoPackage files (*_basinHru.gpkg and *_riverNetwork.gpkg) located in the basin’s shapefiles/ directory. The framework automatically resolves the shapefile directory from the base_model_config_dir path. If GeoPackage files are not found, the map section is silently skipped.

Dependencies: Map generation requires the fiona Python package for reading GeoPackage files:

pip install fiona

Post-Calibration Pipeline

The full post-calibration pipeline runs automatically after DDS optimization completes:

Save metadata: parameter_definitions.json and matched_data.csv are written to results/
Generate per-variant report: calibration_report.html with interactive charts, maps, and parameter tables
Detect sibling workspaces: Scan parent directory for workspace_{basin_id}_* directories
Generate basin report: If multiple completed variants are found, basin_report.html is created in a shared basin_reports/ directory

To enable automatic basin report generation, include basin_id and variant in your calibration configuration:

config = {
    "basin_id": "CAN_01AM001_meso",
    "variant": "A",
    # ... other settings ...
}

The pipeline requires no manual intervention — simply run your calibration and reports appear in the output directory.

Complete Configuration Example

# =============================================================================
# PATH CONFIGURATION
# =============================================================================
base_model_config_dir = "/path/to/Model_Config"
calibration_work_dir = "/path/to/calibration_workspace"
observation_data_path = "/path/to/observations.csv"
test_case_dir = "/path/to/test_case"

# =============================================================================
# CONTAINER RUNTIME
# =============================================================================
container_runtime = "docker"  # or "apptainer"
docker_container_name = "docker_openwq"
docker_compose_path = "/path/to/docker-compose.yml"
executable_name = "mizuroute_lakes_openwq_Release"
executable_args = ""  # No flags - file manager is positional argument
executable_full_path = "/code/.../bin/mizuroute_lakes_openwq_Release"
file_manager_path = "/code/.../mizuroute_in/settings/mizuroute.control"

# mpirun -np 2 required for mizuRoute domain decomposition
command_template = (
    "cd {eval_dir} && mpirun --allow-run-as-root -np 2 "
    "-x master_json {exec_path} {file_manager}"
)

# =============================================================================
# CALIBRATION PARAMETERS
# =============================================================================
calibration_parameters = [
    {
        "name": "k_nitrification",
        "file_type": "bgc_json",
        "path": ["CYCLING_FRAMEWORKS", "N_cycle", "3", "PARAMETER_VALUES", "k"],
        "initial": 0.03,
        "bounds": (0.001, 0.5),
        "transform": "log"
    },
    {
        "name": "fertilizer_N_scale",
        "file_type": "ss_csv_scale",
        "path": {"species": "all"},
        "initial": 1.0,
        "bounds": (0.1, 5.0)
    },
]

# =============================================================================
# CALIBRATION SETTINGS
# =============================================================================
algorithm = "DDS"
max_evaluations = 500
n_parallel = 1
objective_function = "KGE"
objective_weights = {"NO3-N": 1.0, "NH4-N": 0.5}
calibration_targets = {
    "species": ["NO3-N", "NH4-N"],
    "reach_ids": "all",
    "compartments": ["RIVER_NETWORK_REACHES"]
}
random_seed = 42

# =============================================================================
# TEMPORAL RESOLUTION
# =============================================================================
temporal_resolution = "monthly"  # native, daily, weekly, monthly, yearly
aggregation_method = "mean"      # mean, sum, median, min, max

# =============================================================================
# RUN CONTROL
# =============================================================================
run_calibration_flag = True  # Set to False to only validate/prepare files

Dependencies

Required:

numpy
pandas
h5py
matplotlib (for plotting)

Optional:

SALib (for Morris/Sobol sensitivity analysis):
```
pip install SALib
```

Runtime:

Docker (local) or Apptainer/Singularity (HPC)

Troubleshooting

“No observations loaded” error:

Check that your observation file path is correct and the file contains data matching your target species and reach IDs.

“Container not found” error:

Ensure Docker is running (local) or the Apptainer SIF path is correct (HPC).

Calibration converges to bounds:

Consider expanding parameter bounds or using log-transform for rate constants.

HPC jobs fail immediately:

Check SLURM logs in slurm_logs/ for error messages. Verify bind paths and module loads in the job script.

Checkpoint resume fails:

Ensure you’re running from the same working directory. Check that calibration_state.json exists and is not corrupted.