Three calibrations to start from

Three complete configurations, shipped as files rather than as snippets, so that a copy is a copy and not a transcription. Each one loads: a unit test validates all three on every commit, which is what keeps them true after a key is renamed.

Each is an overlay. It carries the search and inherits the catchment, the data and the solver from the project it points at, so the model that is calibrated and the model that is run are one description. The project has to declare the parameters the overlay moves, in [flow] param_list and [flow.param.<id>]; a calibration cannot search over a property the model does not have. Copy one next to your project.toml and run it:

hmp calibrate --check calibration_single_gauge.toml   # nothing solves
hmp calibrate calibration_single_gauge.toml

Pick by what the site actually offers.

Recipe	Use it when	What it decides
One gauge	One discharge record, one property to identify.	The bounds, the metric, and how much of the record is burn-in.
Several targets	A gauge, a piezometer and a lake all constrain the same model.	What share of the cost each target carries.
A published method	The catchment is mapped but poorly gauged.	Almost nothing: the method is named, not retyped.

One parameter, one gauge, one metric

The shape almost every calibration starts from. A hydraulic conductivity is swept over a bounded range and scored against one observed discharge series.

Two lines carry most of the modelling decision. bounds states what the site could plausibly be, and a range far wider than the evidence buys a search over a region the data already excludes. warmup_periods drops the head of the record from the score, because that stretch is still forgetting the initial condition; size it by raising it until the objective stops moving.

# Reference calibration 1 of 3: one parameter, one gauge, one metric.
#
# The shape almost every calibration starts from. One hydraulic property is
# searched over a bounded range, the run is scored against one observed
# discharge series, and the best value is replayed with its full outputs so the
# figures show the calibrated model and not the last trial.
#
# Load it with `hmp calibrate calibration_single_gauge.toml`.

# This file is an overlay: it carries the calibration and inherits the
# catchment, the data and the solver from the project it points at. Every
# calibration in the repository is written this way, so the same model is
# run and calibrated from one description.
base_config = "project.toml"

[workflow]
mode = "calibration"

[simulation.time]
start_datetime = "2000-01-01"
end_datetime = "2005-12-31"
step_value = 1
step_unit = "day"

[calibration]
# 'grid' sweeps the range on a regular lattice: no tuning, reproducible, and it
# shows the shape of the cost rather than a single point. Switch to 'optuna' or
# 'scipy_nelder_mead' once the shape is known and the budget matters.
method = "grid"
max_iter = 25
seed = 42

# Which observed series the cost is computed on, and how.
variable = "discharge"
objective = "nse_log"

# The first year is still forgetting its initial condition. Dropping it from the
# score keeps that memory out of the calibrated value. Size it by raising it
# until the objective stops moving, not by picking a round number.
warmup_periods = 365

# Keep one full simulation: the best trial, replayed with every output.
save_runs = "best_n"
save_best_n = 1
rerun_best_with_outputs = true

[calibration.parameters.K]
# Bounds are the statement of what the site could plausibly be. Six decades is
# what a fractured basement aquifer spans; narrow it when a pumping test or a
# regional map says more.
bounds = [1e-7, 1e-3]
transform = "log"
path = "flow.param.K.field.value"
units = "m/s"

Several targets, weighted against each other

A catchment rarely offers one observation. Here the model is scored at once on the outlet hydrograph, on a piezometer, and on a lake level, each with its own metric and its own share of the cost.

The weights are the decision this file exists to make explicit. They sum to 1.0, so each reads directly as a share: the hydrograph carries 65 %, the piezometer 25 %, the lake 10 %. Write them any way you like; only the ratios matter to the search, and summing to one is what makes them readable.

normalize_cost matters here and did not in the previous recipe. An RMSE on heads is in metres, an NSE cost is dimensionless, and adding them raw lets the unit set the weighting instead of the weight.

observes is what makes the weighting mean anything. It names a station the project loaded; the record is aligned on the simulated timestamps, and the block scores that. Without it a block can only be scored against observed_values, a positional vector transcribed into the file with no dates on it, so the one route that can weight several targets was also the one that could not read a real observation.

Two consequences follow. A scoring_window becomes applicable, because the samples now carry dates to cut on; it is still refused for any output scored positionally, by name. And each block reports <output>.n_paired, how many dated samples the alignment kept: a weight of 65 % resting on three surviving days is not what the file says it is, and that is where it shows.

# Reference calibration 2 of 3: several targets, weighted against each other.
#
# A catchment rarely offers one observation. Here the model is scored at once on
# the outlet hydrograph, on a piezometer, and on a lake level, each with its own
# metric and its own share of the cost. Two parameters move together.
#
# The weights are the modelling decision this file exists to make explicit: they
# say what the calibration is for. Read them as shares of the total cost.

# This file is an overlay: it carries the calibration and inherits the
# catchment, the data and the solver from the project it points at. Every
# calibration in the repository is written this way, so the same model is
# run and calibrated from one description.
base_config = "project.toml"

[workflow]
mode = "calibration"

[simulation.time]
start_datetime = "2000-01-01"
end_datetime = "2005-12-31"
step_value = 1
step_unit = "day"

[calibration]
method = "optuna"
max_iter = 200
seed = 42
# A window in dates rather than warmup_periods in samples: it says the same
# thing whatever the output frequency. The two are mutually exclusive.
# Every scored output names a station, so the samples carry dates and the window
# applies exactly. Without `observes` this would be refused: there would be no
# time axis to cut on. 2000 is the spin-up year, simulated and not scored.
save_runs = "best_n"
save_best_n = 5
rerun_best_with_outputs = true
# 'full' keeps the per-block raw and normalized costs of every trial, which is
# what makes a weighting argument checkable after the fact.
persist_iteration_detail = "full"

[calibration.scoring_window]
start = "2001-01-01"
end = "2005-12-31"

[calibration.parameters.K]
bounds = [1e-7, 1e-3]
transform = "log"
path = "flow.param.K.field.value"
units = "m/s"

[calibration.parameters.Sy]
bounds = [5e-3, 0.35]
transform = "log"
path = "flow.param.Sy.field.value"
units = "-"

# --- What is observed -------------------------------------------------------
#
# `observes` names a station the project loaded. Its record is aligned on the
# simulated timestamps, so a block scores a real chronicle rather than a vector
# transcribed into this file. The data family follows `variable`: discharge from
# hydrometry, head from piezometry, stage from lake_levels.

[calibration.outputs.outlet_discharge]
variable = "discharge"
support = "boundary"
boundary_id = "outlet"
observes = "J7branch"

[calibration.outputs.piezometer_p1]
variable = "head"
support = "point"
x = 348120.0
y = 6812450.0
observes = "BSS001ABCD"

[calibration.outputs.lake_level]
variable = "stage"
support = "lake"
lake_id = "cheze"
observes = "cheze_staff_gauge"

# --- How each target is scored, and for how much -----------------------------
#
# Three blocks summing to 1.0, so each weight reads directly as a share:
# the hydrograph carries 65 %, the piezometer 25 %, the lake 10 %.
#
# normalize_cost matters here and not in the single-target file: an RMSE on
# heads is in metres, an NSE cost is dimensionless, and summing them raw would
# let the unit decide the weighting instead of the weight. It applies to the
# piezometer block alone, because that is the only cost carrying a unit; on an
# efficiency score it is refused rather than silently reweighting the block.
#
# Each block reports `<output>.n_paired`: how many dated samples the alignment
# actually kept. A weight of 65 % on three surviving days is not what the file
# says it is, and that number is where it shows.

[[calibration.objective_blocks]]
name = "hydrograph"
metric = "nse_log"
weight = 0.65
uses_outputs = ["outlet_discharge"]

[[calibration.objective_blocks]]
name = "piezometry"
metric = "rmse"
weight = 0.25
uses_outputs = ["piezometer_p1"]
normalize_cost = true

[[calibration.objective_blocks]]
name = "lake"
# 'reservoir' is half nse_seasonal and half nse_delta: an impounded level is
# beaten by its own seasonal cycle under a plain NSE, and its flux errors show
# in the increments. Both halves are efficiency scores, so there is no unit in
# the cost to remove and normalize_cost is refused here: the share is the weight.
metric = "reservoir"
weight = 0.10
uses_outputs = ["lake_level"]


# ---------------------------------------------------------------------------
# How wide the answer is
# ---------------------------------------------------------------------------
# `multistart` runs the whole search `restarts` times from `restarts` different
# starting points and reports the spread of the optima it reaches. It is the only
# declared method that can see a second basin: `cost_profile` reads one search's
# own trace, so it says what that search could not tell apart and never what it
# never visited.
#
# What it costs is not hidden: four restarts of a two-hundred-evaluation search is
# eight hundred model runs. Start from `cost_profile`, which is free, and pay for
# this one when a parameter matters enough that you need to know whether another
# optimum exists.
#
# The calibrated value does not move: it is the best of the restarts, still one
# number you can write into a report. The spread sits beside it, and a parameter
# whose restarts landed more than a factor ten apart is flagged as not identified
# by this calibration, however tight any single restart looked.
#
# Refused on `grid` and `bisection`: an exhaustive sweep and a root search return
# the same answer every time, so their spread would state a certainty nothing
# established.
[calibration.uncertainty]
method = "multistart"
restarts = 4

A published method, named instead of retyped

Conductivity from the extent of the hydrographic network, storage from the hydrograph. See Calibrating on the stream network for what the criterion measures and [Abhervé et al., 2023] for the method.

protocol = "matching_hydrographic_network" writes the whole assembly: two stages, their criteria, the regime and time-grid overrides that make the first steady and the second transient, and the objective block wiring. What stays in the file is what belongs to the site.

# Reference calibration 3 of 3: a published method, named instead of retyped.
#
# Conductivity from the extent of the hydrographic network, storage from the
# hydrograph. Abherve et al. (2023), doi:10.5194/hess-27-3221-2023.
#
# The whole two-stage assembly is one line. What is left in this file is what
# belongs to the site: the two search ranges and the mapped network. The
# protocol writes the stages, their criteria, and the regime and time-grid
# overrides that make the first steady and the second transient.
#
# `hmp calibrate --list-phases <this file>` prints the assembly before running it.

# This file is an overlay: it carries the calibration and inherits the
# catchment, the data and the solver from the project it points at. Every
# calibration in the repository is written this way, so the same model is
# run and calibrated from one description.
base_config = "project.toml"

[workflow]
mode = "calibration"

[simulation.time]
start_datetime = "1995-01-01"
end_datetime = "2020-12-31"
step_value = 1
step_unit = "day"

# The criterion measures lengths along the flow paths of the DEM, so the mapped
# network has to sit in the talwegs of the routing surface. Burn it first.
[geographic.enforce_streams]
enabled = true
stream_geometry_path = "streams.gpkg"
mode = "constant"
depth_m = 30

[calibration]
seed = 42
save_runs = "best_n"
save_best_n = 1
rerun_best_with_outputs = true
persist_iteration_detail = "full"

# The signed criterion crosses zero once over several decades of K, which is
# what makes the first stage a root search rather than a sweep. The bisection
# refuses a parameter that is not searched in log space.
[calibration.parameters.K]
bounds = [1e-7, 1e-3]
transform = "log"
path = "flow.param.K.field.value"
units = "m/s"

# The lower bound is not the physical floor, deliberately: the solve cost
# climbs steeply as Sy goes to zero, and the observed recession already excludes
# that region.
[calibration.parameters.Sy]
bounds = [5e-3, 0.35]
transform = "log"
path = "flow.param.Sy.field.value"
units = "-"

# The stage-one criterion is solved at zero, so five per cent of the best cost
# is zero and no interval would be reported. State the width in the unit of the
# cost instead: 25 metres of network offset is the agreement this catchment can
# claim, and the interval is the range of K the trials could not tell apart
# inside it.
[calibration.uncertainty]
method = "cost_profile"
mode = "absolute"
tolerance = 25.0

[calibration.outputs.seepage_network]
support = "network"
stream_geometry_path = "streams.gpkg"
# One cell one vote is the paper. Weight by area on a mesh refined along the
# streams, where cell density is highest exactly where distances are smallest.
weighting = "area"
tau_specific_ratio = 1.0e-4
roptim_max = 2.0
# The positional accuracy of the linework does not improve because the mesh
# does. Without this floor, the validity length would follow the mesh.
observed_position_accuracy = "50 m"

# `protocol = "matching_hydrographic_network"` on its own line is the whole
# declaration. It is written here as a table because this file also states the
# options, and TOML will not let one key be both a string and a table.
# Every option below is optional: the defaults reproduce the published method.
[calibration.protocol]
name = "matching_hydrographic_network"
# Which of this file's parameters plays which role.
conductivity = "K"
storage = "Sy"
# Stage one. 'distance_gap' is the signed difference of Eq. 1, whose zero is the
# balance the paper solves for.
steady_metric = "distance_gap"
steady_method = "bisection"
steady_max_iter = 20
# How precisely each stage has to pin its parameter before it stops, as a
# relative precision on that parameter: the paper's one per cent on the
# conductivity, five on the storage. One number per stage, whatever engine walks
# it: each engine's own stopping option is written from it.
steady_tolerance = 0.01
# Stage two. Any registered engine works here: the method is the pair of
# criteria and the order they run in, not the optimizer that walks them.
transient_metric = "nse_log"
transient_method = "scipy_nelder_mead"
transient_max_iter = 120
transient_tolerance = 0.05

# The coarse log sweep the paper runs before bisecting. It checks the
# monotonicity the method assumes instead of supposing it.
[calibration.protocol.steady_optimizer_kwargs]
sweep_points = 7

# The transient stage still has to simulate 1995 to reach 1996 with a state it
# has forgotten its initial condition in. It simulates it and does not score it.
[calibration.protocol.scoring_window]
start = "1996-01-01"
end = "2020-12-31"

A file cannot both name a protocol and declare its own [[calibration.phases]] or [[calibration.objective_blocks]]: that is two answers to one question, and it is refused rather than silently resolved. Drop the protocol to write the stages by hand, or drop the stages to let the protocol write them.

Every option under [calibration.protocol] has a default that reproduces the published method, so the shortest form of this recipe is one line:

[calibration]
protocol = "matching_hydrographic_network"

The engines are not part of the method. steady_method and transient_method take any registered optimizer, so the same two criteria can be walked by a bisection, by Nelder-Mead, or by Optuna without changing what is being calibrated.

Checking before the solver starts

A calibration is hours of solver time, so the last thing you want is a typo found at hour three. One flag runs every static check and solves nothing:

hmp calibrate --check project.toml

ERROR   [calibration.parameters.K]: path 'flow.param.Kh.field.value' is not a
        value this configuration carries. `hmp config targets` lists them.
ERROR   [[calibration.objective_blocks]] 'hydrograph': uses_outputs names ghost,
        which [calibration.outputs] does not declare. Declared: outlet.
calib.toml: 2 error(s), 0 warning(s).

Every check runs, so a file with three mistakes comes back with three findings and takes one pass to fix. It checks that each parameter path is a value the configuration actually carries, that the bounds are ordered and inside the physical range the registry enforces, that a stream geometry is where the run will look for it, and that every name a block or a phase uses is declared. It exits on the config code when anything is wrong, so a script can gate on it.

It also faces each search with what its engine says it can be handed. A bisection moves one parameter, walks a log10 variable and drives a signed residual to zero; hand it two parameters, a linear one, or a search scored on nse_log, and the refusal used to arrive when that phase started, after the phases before it had spent their whole budget.

What it cannot see it does not pretend to: an observed record is loaded by the data step, which needs a delineated catchment, so preflight checks the names a file declares against each other and leaves the loading to the run.

Finding what a project can calibrate

A parameter is declared by a dotted path into the configuration, and guessing one is not a workflow. Ask the project:

hmp config targets project.toml

path                                current      unit  physical range
flow.bc.drainage.value                0.001      m2/s  -
flow.param.K.field.value              5e-05       m/s  1e-14 .. 100
flow.param.Ss.field.value             1e-05       m-1  1e-09 .. 0.001
flow.param.Sy.field.value              0.05         -  0.0001 .. 0.5
flow.sinks_sources.recharge.values         0    mm/day  -

The list comes from the resolved configuration, so it names the parameters and boundaries this project actually declares, not everything the schema could hold. The three columns are what a bound is written from: where the value sits today, in what unit, and the range the physical registry will refuse outside of. A range shown as - means the registry does not know this identifier, so nothing will check the bounds for you.

Add --json for the same catalogue as machine-readable records.

Reading the value the search returns

A calibration returns one number per parameter, and a number on its own says nothing about its own standing. The report adds two things read off the trials the search already ran, so they cost no extra model runs.

How wide the optimum is. parameter_intervals gives, per parameter, the range of sampled values whose cost stayed within 5 % of the best, together with how many trials that was out of how many. Read it as “the search could not tell these apart”, not as a confidence interval: it rests on no error model, and a parameter the search never varied far has a narrow range because nothing else was tried.

The flag to read first is whether the range runs into a search bound. There the record did not determine the parameter, the search simply ran out of room, and widening the bounds is the next step rather than reporting the edge as a result. The run says so in one line:

K = 3.1e-06, and 14 of 60 trials scored within 0.081 of the best over
[2.4e-06, 4.8e-06].
Sy = 0.35, and 31 of 120 trials scored within 0.12 of the best over
[0.19, 0.35]; that range runs into the upper search bound.

The width is a decision, so it is written in the file:

[calibration.uncertainty]
method = "cost_profile"
mode = "relative"
tolerance = 0.05

A criterion whose best cost is zero, such as the stream-network gap, has no fraction of itself to take: five per cent of zero is zero and no interval comes back. State the width in the unit of the cost there instead, and the run says so when it has to:

[calibration.uncertainty]
mode = "absolute"
tolerance = 25.0        # metres of network offset

How the members were made addable. [calibration.aggregate] names it rather than leaving it implicit:

[calibration.aggregate]
weighting = "manual"        # or "error": one over sigma
nested_gauges = "total"     # or "incremental" on imbricated catchments
min_samples = 30            # refuse a member scored on fewer pairs
on_member_failure = "veto"  # or "drop", which records what was left out

Two questions the word “weight” runs together. Whether an error is large for what the instrument can resolve is a property of the measurement, not a decision. What matters more between the outlet and the reservoir is a decision, and yours. The cost is the product of both.

weighting = "error" divides each residual by what its gauge resolves, so the members become pure numbers before the weights apply. It needs a residual metric and a loaded record carrying an error model, and is refused without both: an efficiency score has no residual to divide, and a vector typed into the file carries no sigma.

nested_gauges matters when two gauges sit on imbricated catchments. Nothing is double counted, but the residuals are statistically dependent and no standard correction exists. "total" scores each gauge against its own full drained area, which is what a gauge measures; "incremental" scores the downstream one on what its own reach adds, which is the only mechanisable way to make the two independent. The overlap is measured and reported either way.

Whether two parameters were told apart. correlated_parameters lists the pairs that moved together across the whole search to hold the same cost. Such a pair was not identified: the search stopped somewhere on a ridge and reported that point as a minimum. The sign is kept, because it says which way the trade-off ran.

Why the mesh is not a parameter

A mesh resolution or a refinement setting cannot be declared under [calibration.parameters]; the file is refused with the reason. The stream-network criterion is normalised by cell size, so refining the mesh moves the yardstick the search is scored against, and a search that optimises it improves the number by changing the ruler.

A mesh question is a convergence question, and it is answered by a sweep: one run per mesh, compared, and read for convergence rather than for a best score.

# mesh_sweep.toml
[workflow]
mode = "comparison"

[comparison]
comparison_id = "mesh_convergence"
base_simulation_config = "project.toml"
output_root = "outputs/mesh_convergence"
reference_simulation = "mesh_250"

[[comparison.simulation]]
id = "mesh_500"
label = "target cell size 500 m"
solver = "modflow6"

[comparison.simulation.overlay.mesh_catchment.zone_meshing]
global_size = 500.0

[[comparison.simulation]]
id = "mesh_350"
label = "target cell size 350 m"
solver = "modflow6"

[comparison.simulation.overlay.mesh_catchment.zone_meshing]
global_size = 350.0

[[comparison.simulation]]
id = "mesh_250"
label = "target cell size 250 m"
solver = "modflow6"

[comparison.simulation.overlay.mesh_catchment.zone_meshing]
global_size = 250.0

# What the three meshes are compared on. The seepage extent is the quantity
# the stream-network criterion reads, so it is the one whose mesh sensitivity
# decides how fine the calibration has to be.
[[comparison.observable]]
name = "seepage_map_last"
variable = "seepage_areas"
support = "map"
time = "last"
unit = "-"

[[comparison.observable]]
name = "head_map_last"
variable = "watertable_elevation"
support = "map"
time = "last"
unit = "m"

hmp compare mesh_sweep.toml

A cell may be a whole calibration rather than a single run, which is what a structural sweep actually is: one complete calibration per mesh, each with its own plan. Declare it in the cell’s overlay:

[comparison.simulation.overlay.workflow]
mode = "calibration"

[comparison.simulation.overlay.calibration]
protocol = "matching_hydrographic_network"

The child is materialised the same way and dispatched by hmp run on its own [workflow] mode, so the cohort machinery does not need to know which kind of cell it built. What it must not do is rank them: the criterion is normalised by cell size, so the cheapest cost belongs to the finest mesh whatever the hydrogeology.

Read the spread across the three as the numerical error on whatever you report, and refine until it stops moving. Calibrate on the coarsest mesh whose answer no longer changes: a finer one buys solve time, not information. Running the calibration on each mesh in turn and keeping the K that scored best is the same circularity in a longer form, because the criterion those K are compared on is itself normalised by cell size.