River Flow and Nutrient Load Projection Guide

river-projection (stage 1, below) lives in ocean-prep, alongside bc-correct/bc-diagnostics — install and run it from the ocean-prep repo (cli/river_projection.py), with its config at ocean-prep/config/river_projection_example.yaml. river-diagnostics (stage 2) stays in ocean-post. Confirmed 2026-10-07; this page previously referenced a config/river_projection.yaml in ocean-post that was never there, and the real tool and config were found in ocean-prep instead. No project-specific copy of the config exists yet in this checkout — the only one found is ocean-prep’s own example.

Overview

This guide describes how to project future river flows and nutrient loadings using the delta-change method: bias-corrected gridded precipitation (P) and evapotranspiration (E) from CMIP6 are used to scale the observed seasonal discharge cycle from the EMORID dataset forward in time. Nutrient loads are then computed by multiplying projected discharge by an observed flow-weighted mean concentration.

EMORID observed Q + nutrients
        │
        ├─── monthly Q climatology ──────────────────────────┐
        │    (calibration period)                            │
        │                                                    ▼
ERA5 P, E (reference)                            Q_future = Q_clim(month) × ΔPE(t)
        │                                                    │
        └─── P–E climatology ────────────────► ΔPE(t)       │
             (calibration period, denominator)               │
                                                             │
BC future P, E (CMIP6)                                       ▼
        │                                        nutrient loads = Q_future × FWMC
        └─── future P–E ──────────────────────► (numerator)

Commands

# Stage 1 (ocean-prep): project future flows and nutrient loads
river-projection --config config/river_projection_example.yaml --model MPI-ESM1-2-HR --scenario ssp245
river-projection --config config/river_projection_example.yaml --model MPI-ESM1-2-HR --scenario ssp585
river-projection --config config/river_projection_example.yaml --model MPI-ESM1-2-HR --scenario ssp245 --dryrun
river-projection --config config/river_projection_example.yaml --model MPI-ESM1-2-HR --scenario ssp245 --allow-missing

# Stage 2 (ocean-post): generate diagnostic plots and tables from the projection output
river-diagnostics --config config/river_projection_example.yaml --model MPI-ESM1-2-HR --scenario ssp245
river-diagnostics --config config/river_projection_example.yaml --model MPI-ESM1-2-HR --scenario ssp245 --analyses-dir /data/local/analyses

--model is required by river-projection (either metadata.model in the config or --model) — it resolves {model} in the forcing and output paths alongside --scenario.

Config on GitHub: config/river_projection_example.yaml (in ocean-prep)


Input data

DatasetLocationContent
EMORID observed flows${EMORID_FOLDER}/to_netcdf/EMORID_1990_2024.ncdaily Q + nutrients
ERA5 precipitation${ERA5_FOLDER}/era5_tp_{year}.ncHourly accumulation (m), variable tp — reference P–E climatology
ERA5 evaporation${ERA5_FOLDER}/era5_e_{year}.ncHourly accumulation (m), variable e, negative for upward flux
BC future P/E (separate)from bc-correct (ocean-prep)Daily, one file per year, variables pr/evspsbl (kg m⁻² s⁻¹)
BC future P−E (composite)from bc-correct (ocean-prep)Jointly corrected pe, preferred over separate pr/evspsbl when present — see “The pe composite” below

EMORID variables available: Q, TotalN, NO3, NH4, DIN, TotalP, PO4, Si, DOC, DIC, SPM. Countries covered: 18 (NE Atlantic and Baltic rim).


Method

Step 1 — Station filtering

Stations are filtered to retain only those with ≥ min_years (default 10) of valid daily Q. This removes gauges with very short or mostly-missing records that would produce unreliable climatologies.

Stations in Norway, Finland, Sweden, and Iceland are flagged with regulated_flag = True in the output but are not excluded. Rivers in these countries are heavily influenced by snowmelt timing and hydropower operations, so the P–E relationship does not hold as reliably as for rainfall-dominated catchments. Projections for flagged stations should be interpreted with caution.

Step 2 — Associate stations with gridded P–E

Each station is paired with the nearest grid cell in the ERA5 P and E fields (for the climatology) and the bias-corrected future P and E fields (for the projection) using nearest-neighbour selection on (lat, lon).

This is a practical approximation: in the absence of catchment boundary polygons the point extraction captures the local climate signal but not the integrated catchment response. A future improvement would ingest catchment boundary polygons (e.g. HydroSHEDS) and compute proper area-weighted P–E over each drainage basin. The delta-change method is robust to this approximation because it relies only on the relative change in P–E, not its absolute value.

Step 3 — Monthly P–E climatology (calibration period)

The mean monthly P–E at each station is computed over the calibration period (e.g. 1990–2019) from ERA5, which is the observational reference:

PE_ref_clim(m) = mean over t in calibration period of [P_ERA5(t) - E_ERA5(t)]
                 where month(t) = m

This 12-value climatology forms the denominator of the delta-change factor.

Why ERA5 and not bias-corrected historical? The BC future files are anchored to the ERA5 climatological distribution by construction. Dividing the future P–E by the ERA5 climatology therefore recovers the pure CMIP6 climate-change signal. Bias-corrected historical files would give the same result (ERA5 ≈ BC-historical by construction) and are not needed.

ERA5 units and sign convention: ERA5 tp/e are hourly accumulations in metres; CMIP6’s bias-corrected pr/evspsbl are daily-mean fluxes in kg m⁻² s⁻¹. Both ERA5 variables are resampled to daily totals (resample: daily, summing the 24 hourly values) and then scaled by 0.011574 (= 1000/86400/24) to match that flux convention. ERA5 e is also stored as a negative accumulation (upward flux = surface loss), so its scale factor is -0.011574 — the sign flip and the unit conversion combined into one number.

Step 4 — Observed monthly Q climatology (calibration period)

The mean monthly observed discharge at each station is computed from EMORID over the same calibration period:

Q_clim(m) = mean over t in calibration period of Q_obs(t)
             where month(t) = m

This is the baseline seasonal cycle that the delta-change method preserves.

Step 5 — Delta-change factor

For each future time step t the delta-change factor is:

ΔPE(t) = PE_future(t) / PE_hist_clim(month(t))

where PE_future(t) = P_future(t) − E_future(t) from the bias-corrected future forcing.

Interpretation: ΔPE(t) = 1.2 means the future P–E in that month is 20 % higher than the historical climatological value, so the projected discharge is also 20 % higher.

Handling non-positive P–E: In arid regions or summer months, the historical P–E climatology can be zero or negative (evaporation exceeds precipitation). Division by a non-positive value is physically undefined for this method; those cells are set to ΔPE = 1.0 (no change) and logged as warnings. If many stations trigger this condition, consider extending the calibration period or reviewing the evaporation variable.

Step 6 — Projected river discharge

Q_future(t) = Q_clim(month(t)) × ΔPE(t)

Values are clipped to zero (no negative discharge).

Why delta-change rather than a full rainfall-runoff model? The delta-change method preserves the observed seasonal cycle and inter-annual variability pattern of discharge without requiring catchment area, soil parameters, or routing. It applies only the climate-driven fractional change in water availability (P–E) to the observed flow. The method is widely used for first-order assessments when catchment characteristics are not available.

Step 7 — Flow-weighted mean concentration (FWMC)

For each nutrient species and station:

FWMC = Σ(Q(t) × C(t)) / Σ(Q(t))   over the FWMC period

The FWMC period is configured independently from the calibration period. Use a recent sub-period (e.g. 2010–2023) rather than the full record to reflect current nutrient conditions; older data from the 1990s may overestimate future loads in rivers that have improved under EU nutrient directives such as the Water Framework Directive and Nitrates Directive.

Both Q and C must be finite for a time step to enter the sum. Stations with no valid Q–C overlap for a given nutrient receive FWMC = NaN and are excluded from load projections for that nutrient.

Step 8 — Projected nutrient loads

Load_future(t) = Q_future(t) × FWMC   [kg s⁻¹ if C in kg m⁻³]

This assumes in-catchment nutrient sources and in-stream processing remain at their current level. This is a conservative baseline; a separate scenario analysis can relax it by applying trend-based or policy-driven adjustments to FWMC.


Configuration

# ocean-prep/config/river_projection_example.yaml

metadata:
  scenario: ssp126           # override at run time with --scenario
  model:    MPI-ESM1-2-HR    # override at run time with --model

emorid:
  path: "${EMORID_FOLDER}/to_netcdf/EMORID_1990_2024.nc"

station_filter:
  min_years: 10        # minimum years of valid daily Q

calibration:
  start: "1990-01-01"
  end:   "2019-12-31"

future:
  start: "2015-01-01"
  end:   "2099-12-31"

fwmc:                  # period for flow-weighted mean concentrations
  start: "2010-01-01"  # omit this block to use the full EMORID record
  end:   "2023-12-31"

nutrients:
  - TotalN
  - TotalP
  - DIN

domain:
  lon_min: -20.0
  lon_max:  40.0
  lat_min:  45.0
  lat_max:  75.0

forcing:
  # ERA5: denominator of the delta-change factor (P–E climatology).
  reference:
    pr:
      path:         "${ERA5_FOLDER}/era5_tp_{year}.nc"
      variable:     tp
      scale_factor: 0.011574    # hourly accumulation (m) → daily kg m⁻² s⁻¹
      resample:     daily       # hourly → daily sum
    evspsbl:
      path:         "${ERA5_FOLDER}/era5_e_{year}.nc"
      variable:     e
      scale_factor: -0.011574   # sign flip (ERA5 e is negative) + unit conversion
      resample:     daily

  # BC future: numerator of the delta-change factor. {model}/{scenario}
  # substituted from metadata/--model/--scenario; {method} (regridding
  # method, e.g. qdm/bilinear) is left as a glob (*) rather than hardcoded,
  # since pr and evspsbl may use different methods.
  future:
    pr:
      path:     "${BIAS_CORRECTED_FOLDER}/CMIP6/{model}/{scenario}/meteo/pr_bc_*_{scenario}_{year}.nc"
      variable: pr
    evspsbl:
      path:     "${BIAS_CORRECTED_FOLDER}/CMIP6/{model}/{scenario}/meteo/evspsbl_bc_*_{scenario}_{year}.nc"
      variable: evspsbl
    # Jointly bias-corrected P-E composite, preferred over separate pr −
    # evspsbl when present (auto-detected by real file presence per
    # model/scenario — not every model produces it, e.g. GFDL-ESM4's CMIP6
    # hfls is Amon-only). See "The pe composite" below.
    pe:
      path:     "${BIAS_CORRECTED_FOLDER}/CMIP6/{model}/{scenario}/meteo/pe_bc_*_{scenario}_{year}.nc"
      variable: pe

output:
  dir:          "${BIAS_CORRECTED_FOLDER}/CMIP6/{model}/{scenario}/rivers"
  compress:     true
  analyses_dir: "${OCEANICU_ANALYSES_FOLDER}"  # not used by river-projection itself today —
                                                 # a placeholder for river-diagnostics (stage 2)

{model} and {scenario} in the future forcing paths and output.dir are substituted at runtime from metadata.model/metadata.scenario (or --model/--scenario). This makes it straightforward to run multiple models and SSPs from the one shared config:

for scenario in ssp126 ssp370 ssp585; do
    river-projection --config config/river_projection_example.yaml \
        --model MPI-ESM1-2-HR --scenario $scenario
done

The pe composite

When forcing.future.pe is present in the config, river-projection probes for real output at that path for the current --model/--scenario before using it: if found, it uses the jointly bias-corrected P−E composite instead of separately-corrected pr − evspsbl (separate correction does not guarantee the difference preserves ERA5’s own P–E climatology). If no pe output exists for that model/scenario, it logs this and falls back to pr/evspsbl automatically — the same shared config works for every model via --model without a per-model variant.


Output

Two NetCDF files are written per run, both with dimensions (site, time):

FileVariablesDescription
river_flows_future_<scenario>.ncQ, regulated_flagProjected monthly discharge (m³ s⁻¹)
nutrient_loads_future_<scenario>.ncTotalN, TotalP, DIN, …Projected loads (kg s⁻¹)

Global attributes record the scenario name, calibration and future periods, and the FWMC period. These are written to output.dir (not into the analyses_dir/validations tree) — output.analyses_dir in the config is not used by river-projection itself; it only matters to river-diagnostics below.


Diagnostics

After running river-projection, use river-diagnostics (ocean-post) to generate plots and tables from the projection output. It reads the same config file (for metadata, emorid, and output.dir) and writes into the canonical analyses layout:

<analyses_dir>/scenarios/<model>/<scenario>/river/
    plots/
        q_change_map_<scenario>.png      # station map coloured by % Q change
        q_time_series_<scenario>.png     # annual mean Q: calibration + near/mid/far future
        q_seasonal_<scenario>.png        # monthly Q climatology by period
        nutrient_loads_<scenario>.png    # annual nutrient load time series
    tables/
        river_projection_statistics_<scenario>.csv
        river_projection_statistics_<scenario>.txt
river-diagnostics --config config/river_projection_example.yaml --model MPI-ESM1-2-HR --scenario ssp245
river-diagnostics --config config/river_projection_example.yaml --model MPI-ESM1-2-HR --scenario ssp585
river-diagnostics --config config/river_projection_example.yaml --model MPI-ESM1-2-HR --scenario ssp245 --analyses-dir /data/local/analyses

--analyses-dir overrides output.analyses_dir from the config and has no ./analyses fallback — see Tidal Analysis: analyses_dir has no default. --scenario and --model override the corresponding metadata keys.

The diagnostics are automatically picked up by ocean-reporting and included in the Hugo scenarios page under a River Flow and Nutrient Load Projections section.


Dry run

Use --dryrun to check the config and probe forcing-file existence (every probed year, for both ERA5 reference and BC future — plus the pe auto- detect probe when configured) without loading any data:

river-projection --config config/river_projection_example.yaml --model MPI-ESM1-2-HR --scenario ssp245 --dryrun

Both --dryrun and a real run exit non-zero if any probed year is missing — a config/filename mismatch is meant to fail loudly, not silently degrade the projection. Pass --allow-missing only for genuinely partial, expected coverage (e.g. a scenario whose later years are not produced yet).


Limitations and future improvements

LimitationImpactPotential improvement
Nearest-neighbour P–E extractionIgnores catchment shape; grid noise can affect individual stationsIngest HydroSHEDS catchment polygons and compute area-weighted P–E
Constant FWMCDoes not capture land-use change or emission trendsFit a trend to historical C(t) and extrapolate; or apply emission scenario adjustments
Regulated rivers flagged but not excludedDelta-change projections less physically meaningful for snowmelt / hydropower dominated basinsSeparate treatment using temperature-based snowmelt model
No uncertainty from FWMCSingle FWMC value per stationBootstrap confidence intervals on FWMC from the observed record

Key references

  • Arnell (1999), Climate change and global water resources, Global Environmental Change 9, S31–S49. — Foundation of the delta-change approach for river flow projections.
  • Middelkoop et al. (2001), Impact of climate change on hydrological regimes and water resources management in the Rhine basin, Climatic Change 49, 105–128. — Applied delta-change example.
  • Kronvang et al. (2009), Nitrogen and phosphorus losses from agricultural areas in European river basins, Science of the Total Environment. — FWMC methodology context for riverine nutrient loads.