Model Validation

---

Lecture 12

October 7, 2026

Overview

• Review of Last Class
• When Should We Trust A Model?
• Model Calibration: Fitting to Data
• Positive and Negative Controls
• Validation: Trusting A Model
• Key Takeaways and Upcoming Schedule

Review of Last Class

Gaussian Plume Behavior Depends On A Number of Factors

The ground-level peak is downwind from the stack because the plume has to spread down to the ground first.

A taller stack lowers the peak and pushes it further out: the same emissions, spread over more distance.

Atmospheric stability can change the peak, the distance from the stack, and the ground-level area in exceedance of a standard.

Questions?

Text: VSRIKRISH to 22333

URL: https://pollev.com/vsrikrish

See Results

When Should We Trust A Model?

Sources of Model Error

A model will always give you a number. Someone asks: why should I believe this?

1. Structural: the model is the wrong model for this system or question.
2. Parametric: the structure is right, the parameter values are wrong.
3. Numerical: the equations are right, error comes from discretization, Monte Carlo, etc.

Model Calibration: Fitting to Data

Calibrating the DO Model

The rate constants \(k_a\), \(k_c\), and \(k_n\) are approximations, not the “true” decay or reaeration dynamics.

In practice we need to have some data to estimate these parameters in the context of our model. Measuring these rates in a lab may be helpful, or it might not be, depending on how simplified or aggregated the modeled process is.

Example Data We Might Collect

# Seeded so the lecture reports the same numbers every time it renders.
Random.seed!(3)

observation_step = 8                      # one observation every 2 km on a 0.25 km grid
truth = simulate_do(C0, river_length, Δx, true_params)
observation_index = 1:observation_step:length(truth)
observation_distance = [(i - 1) * Δx for i in observation_index]
observations = truth[observation_index] + rand(Normal(0, 0.2), length(observation_index))

p_data = scatter(observation_distance, observations, color=cb_vermillion, markersize=7,
    xlabel="Distance downstream  [km]", ylabel="DO  [mg/L]",
    label="measurements", legend=:bottomright)
plot!(p_data, 0:Δx:river_length, truth, color=:black, linestyle=:dash, linewidth=4, label="'true' model")
plot!(p_data, size=(1100, 430), left_margin=12mm, right_margin=8mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 1: Twenty-one dissolved oxygen measurements, every 2 km, against the model that generated them.

Scoring Rules

rmse(simulated, observed) = sqrt(mean((simulated .- observed).^2))

function do_rmse(ka, kc, kn)
    candidate = (; ka, kc, kn, B0=true_params.B0, N0=true_params.N0,
                   Cs=true_params.Cs, U=true_params.U)
    simulated = simulate_do(C0, river_length, Δx, candidate)
    return rmse(simulated[observation_index], observations)
end

A scoring rule scores how closely the model matches the measurements. A common (and typically good!) choice:

\[\text{RMSE} = \sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(y^\text{sim}_i - y^\text{obs}_i\right)^2}\]

Lower is better; zero is impossible in real life. The true rates score 0.18; the model with \(k_a=k_n=k_b=0.5\) scores 1.02.

Optimizing RMSE

lower_bounds = [0.05, 0.05, 0.05]
upper_bounds = [2.0, 2.0, 2.0]
initial_guess = [0.5, 0.5, 0.5]

calibration = Optim.optimize(v -> do_rmse(v[1], v[2], v[3]),
    lower_bounds, upper_bounds, initial_guess)
fitted = Optim.minimizer(calibration)
@show fitted;

fitted = [0.6106236141326572, 0.34699675218765863, 0.3469967525444313]

The Calibrated Fit

fitted_params = (; ka=fitted[1], kc=fitted[2], kn=fitted[3],
                   B0=true_params.B0, N0=true_params.N0,
                   Cs=true_params.Cs, U=true_params.U)

p_fit = scatter(observation_distance, observations, color=cb_vermillion, markersize=7,
    xlabel="Distance downstream  [km]", ylabel="DO  [mg/L]",
    label="measurements", legend=:bottomright)
plot!(p_fit, 0:Δx:river_length, truth, color=:black, linestyle=:dash, linewidth=4, label="true model")
plot!(p_fit, 0:Δx:river_length, simulate_do(C0, river_length, Δx, fitted_params),
    color=cb_blue, linewidth=3, label="calibrated model")
plot!(p_fit, size=(1100, 430), left_margin=12mm, right_margin=8mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 2: The calibrated model against the measurements. It passes through the data at least as well as the truth does.

We Didn’t Recover the “True” Rates

	True	Calibrated
\(k_a\) reaeration	\(0.60\)	0.61
\(k_c\) CBOD decay	\(0.40\)	0.35
\(k_n\) NBOD decay	\(0.25\)	0.35
RMSE (mg/L)	0.18	0.16

And this is when the data-generating model had our model’s form: not the case in reality.

Why Couldn’t We Recover The “Truth”?

cbod_rates = 0.10:0.01:0.80
nbod_rates = 0.05:0.01:0.60
error_surface = [do_rmse(0.6, kc, kn) for kn in nbod_rates, kc in cbod_rates]

p_valley = contourf(cbod_rates, nbod_rates, error_surface, color=:viridis, linewidth=0,
    colorbar_title="\nRMSE  [mg/L]", colorbar_titlefontsize=18,
    xlabel="CBOD decay rate  [1/d]", ylabel="NBOD rate  [1/d]")
scatter!(p_valley, [0.4], [0.25], color=:white, markersize=11, markershape=:star5,
    label="truth", legend=:topright)
scatter!(p_valley, [fitted[2]], [fitted[3]], color=cb_vermillion, markersize=9,
    label="calibrated")
plot!(p_valley, size=(1200, 300), left_margin=12mm, right_margin=22mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 3: RMSE in mg/L over the two decay rates, with reaeration held at its true value. The minimum is a trough, not a point.

The data only constrain the sum of two decays. These two parameters are non-identifiable.

How Does Noise Influence the Fit?

# Score a candidate against any set of measurements, not just the ones above.
function rmse_against(ka, kc, kn, measured)
    candidate = (; ka, kc, kn, B0=true_params.B0, N0=true_params.N0,
                   Cs=true_params.Cs, U=true_params.U)
    simulated = simulate_do(C0, river_length, Δx, candidate)
    return rmse(simulated[observation_index], measured)
end

# Nothing about the river changes across these refits. The only difference is
# which measurements we happened to take.
n_refits = 8
refit_ka = zeros(n_refits)
refit_kc = zeros(n_refits)
refit_kn = zeros(n_refits)

for seed in 1:n_refits
    Random.seed!(seed)
    draw = truth[observation_index] + rand(Normal(0, 0.2), length(observation_index))
    result = Optim.optimize(v -> rmse_against(v[1], v[2], v[3], draw),
        lower_bounds, upper_bounds, initial_guess)
    best = Optim.minimizer(result)
    refit_ka[seed] = best[1]
    refit_kc[seed] = best[2]
    refit_kn[seed] = best[3]
end

Rate	True	Lowest	Highest	Spread
\(k_a\) reaeration	\(0.60\)	0.58	0.61	0.03
\(k_c\) CBOD decay	\(0.40\)	0.32	0.45	0.13
\(k_n\) NBOD decay	\(0.25\)	0.2	0.35	0.15

Refitting to new noise moves the decay rates along the trough: the data constrain a combination of the rates.

Why Compensating Processes Make a Valley

Both oxygen sinks decay exponentially in distance:

\[k_c B_0 e^{-k_c x/U} + k_n N_0 e^{-k_n x/U}\]

Over the measured stretch, raising \(k_c\) and lowering \(k_n\) leaves this sum nearly unchanged: the data constrain the sum, not the terms. Reaeration pulls toward \(C_s\) instead, so \(k_a\) stays identifiable.

Positive and Negative Controls

We Just Ran a Positive Control

Data from known parameters: a test with a known right answer, run on the procedure rather than on the system.

• The procedure fits the data. It does not recover the parameters.
• The failure is not a bug. It is a property of this model and this data.
• We only know that because we knew the answer in advance. Real data is always messier.

Negative Control: Fitting Pure Noise

# Same procedure, but the "measurements" are noise about a flat line ---
# there is no sag here for any parameter set to explain.
Random.seed!(11)
flat_observations = fill(mean(observations), length(observation_index)) +
    rand(Normal(0, 0.2), length(observation_index))

function noise_rmse(ka, kc, kn)
    candidate = (; ka, kc, kn, B0=true_params.B0, N0=true_params.N0,
                   Cs=true_params.Cs, U=true_params.U)
    rmse(simulate_do(C0, river_length, Δx, candidate)[observation_index], flat_observations)
end

noise_fit = Optim.minimizer(Optim.optimize(v -> noise_rmse(v[1], v[2], v[3]),
    lower_bounds, upper_bounds, initial_guess))
@show round.(noise_fit; digits=2);

round.(noise_fit; digits = 2) = [0.41, 0.1, 0.19]

A negative control feeds the procedure something with no signal in it. Handed noise, it still returns a single best-fitting set of rates, just as it did for the DO measurements.

An optimiser always returns its argument minimum. Getting an answer is not evidence that there was one.

The Model Has A Structural Sag

noise_params = (; ka=noise_fit[1], kc=noise_fit[2], kn=noise_fit[3],
                  B0=true_params.B0, N0=true_params.N0,
                  Cs=true_params.Cs, U=true_params.U)
noise_curve = simulate_do(C0, river_length, Δx, noise_params)
noise_fit_rmse = noise_rmse(noise_fit[1], noise_fit[2], noise_fit[3])
flat_rmse = rmse(fill(mean(flat_observations), length(flat_observations)), flat_observations)

p_noise = scatter(observation_distance, flat_observations, color=cb_vermillion, markersize=7,
    xlabel="Distance downstream  [km]", ylabel="DO  [mg/L]",
    label="noise 'measurements'", legend=:topright)
plot!(p_noise, 0:Δx:river_length, noise_curve, color=cb_blue, linewidth=3,
    label="model calibrated to noise")
plot!(p_noise, size=(1100, 430), left_margin=12mm, right_margin=8mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 4: The model calibrated to pure noise, against the noise it was fitted to.

The sag comes from the model, not the data: DO enters at 6.2 mg/L and the decay terms pull it down before reaeration recovers it, whatever the rates. The fit is also poor: RMSE 0.73 mg/L, against 0.23 for a flat line through the noise and 0.16 for the real calibration. Looking at the fit, not only at the rates the optimiser returned, is what exposes it. Checked on five other noise draws: the sag was 2.4 to 2.5 mg/L deep and the RMSE 0.58 to 0.72 mg/L every time.

Validation: Trusting A Model

Calibration Is Not Validation

• Calibration: choosing parameters to match observations.
• Validation: showing the model is a legitimate basis for your claim (Oreskes et al., 1994).
  • Empirical validity: does it predict data it has never seen?
  • Face validity: are the processes and parameters plausible?

These can fail independently: a model can predict well for the wrong reasons, or be partly right and still predict badly.

A Plume With a Missing Sink

Our Gaussian plume model had no sinks: everything emitted stays in the air.

Real SO\(_2\) is removed in transit: oxidation to sulfate (about 2% per hour) and deposition to dew-wet ground (fast under a shallow night-time plume).

The model is the right kind of model, a Gaussian plume on a stable night, but it leaves out a process: nothing removes SO2 between the stack and the receptor. In reality SO2 is oxidised to sulfate, slowly, at a few percent per hour, and it deposits onto the ground. Deposition is the bigger sink here: a deposition velocity of about a centimetre per second into a plume only about fifty metres deep removes it quickly. Deposition at night can be lower when plants close their stomata, and higher onto dew-wet surfaces, so treat the exact rate as an assumption for the example. The point is that this is a structural error: no value of the model’s parameters adds a sink.

Can We Calibrate Our Way Out?

function wind_rmse(wind_speed)
    trial = (; plume_params..., wind=wind_speed)
    modelled = [ground_concentration(monitor_downwind[i], monitor_crosswind[i], class_F, trial)
        for i in 1:n_monitors]
    return rmse(modelled, monitor_observed)
end

wind_fit = Optim.optimize(v -> wind_rmse(first(v)), [0.3], Float64[20], [2.5])
fitted_wind = first(Optim.minimizer(wind_fit))

Wind speed is uncertain, and the whole field scales as \(1/u\): the obvious knob. Calibrated wind: 3.2 m/s. RMSE falls from 287 to 115 µg/m³ (38% to 15% of the mean observation).

That’s a decent fit. But the true wind speed was 2.5 m/s. The wind speed calibration is compensating for a mixing component in the model structure.

The Calibrated Plume

calibrated_plume = (; plume_params..., wind=fitted_wind)
centreline_km = 0.3:0.05:7.5
true_centreline = [ground_concentration(x, 0, class_F, plume_params; loss=so2_loss)
    for x in centreline_km]
calibrated_centreline = [ground_concentration(x, 0, class_F, calibrated_plume)
    for x in centreline_km]
on_centreline = monitor_crosswind .== 0       # the other monitors sit off to the side

p_plume_fit = scatter(monitor_downwind[on_centreline], monitor_observed[on_centreline],
    color=cb_vermillion, markersize=7,
    xlabel="Distance downwind  [km]", ylabel="Concentration  [µg/m³]",
    label="monitors on the centreline", legend=:topright, ylims=(0, 2000))
plot!(p_plume_fit, centreline_km, true_centreline, color=:black, linestyle=:dash,
    linewidth=4, label="true model")
plot!(p_plume_fit, centreline_km, calibrated_centreline, color=cb_blue, linewidth=3,
    label="calibrated model")
plot!(p_plume_fit, size=(1100, 430), left_margin=12mm, right_margin=8mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 5: The calibrated plume against the monitors on the plume centreline.

Only the 6 centreline monitors are plotted: the model’s curve is different at each crosswind offset, so the other 6 would not sit on it. At plot scale the calibrated plume looks close. Look at the two ends: it is below the true plume near the stack (977 against 1123 µg/m³ at 1.5 km) and above it far out (748 against 547 at 7 km). The residuals on the next slide show that tilt directly.

The Residuals Reveal The Problem

fitted_params = (; plume_params..., wind=fitted_wind)
fitted_modelled = [ground_concentration(monitor_downwind[i], monitor_crosswind[i], class_F,
        fitted_params) for i in 1:n_monitors]
fitted_error = 100 .* (fitted_modelled .- monitor_observed) ./ monitor_observed   # [%]

p_resid = scatter(monitor_downwind, fitted_error, color=cb_vermillion, markersize=8,
    xlabel="Distance downwind  [km]", ylabel="Residual  [%]")
hline!(p_resid, [0], color=:black, linewidth=2, linestyle=:dash)
plot!(p_resid, size=(1100, 300), left_margin=12mm, right_margin=8mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 6: Residuals of the calibrated model (model minus observation, as a percentage of the observation) at the twelve monitors it was fitted to.

RMSE: 15% of the mean. But these aren’t random, there’s a pattern in the tilt: 17% low at 0.5 km, 40% high at 7 km.

If the model were right and only the data were noisy, these points would scatter by about eight percent either side of zero. Instead they tilt: low near the stack, high far out. RMSE squares and averages the errors, which throws away their sign and their order, so it cannot show this. A tilt that grows with distance points to a process that grows with travel time, which here is the missing loss. No wind speed can fix it, because changing the wind scales every point by the same factor. The plot tells us the model is wrong, but not how badly it predicts where we would use it; the next slide asks the question that matters for using it, whether it predicts data it has not seen.

External Validity: Out Of Sample Data

Monitor	Observed (µg/m³)	Calibrated model (µg/m³)	Error
8 km	531	663	+25%
10 km	270	537	+99%

The inflated wind dilutes the whole plume by one constant factor; the sinks remove more SO\(_2\) the further it travels.

Structural errors can often be revealed by looking at residuals or data withheld from the calibration.

So What Would Have Caught It?

Calibration made the fit good (as measured by the score) and left the model wrong.

• Face check: Are all the SO\(_2\) sources and sinks in the model? Do we care for the problem?
• Residual check: a pattern in the fit’s residuals usually means something is missing; one summary number cannot show it.
• Empirical check: test the calibrated model on data it was not fitted to. Only this tells you whether its predictions hold.

Key Takeaways and Upcoming Schedule

Key Takeaways

• A model that fits better can be more wrong. Calibration minimises a metric, not parameter error.
• Nonuniqueness is normal: data often constrain only a combination of parameters.
• A good fit is not validation. Check calibrated parameters against what you know, and test on data you did not fit.
• Controls test the procedure, not the system.
• Empirical and face validity fail independently.

Next Classes

Monday: Fall Break — no class.

Wednesday, Oct 14: Decision models and linear programming.

Assessments

HW5: assigned Wednesday Oct 14, after Fall Break, due Thursday Oct 22.

Quiz 3: Wednesday Oct 21.

References

Oreskes, N., Shrader-Frechette, K., & Belitz, K. (1994). Verification, validation, and confirmation of numerical models in the Earth sciences. Science, 263, 641–646. https://doi.org/10.1126/science.263.5147.641