More Dissolved Oxygen and Simulation

Lecture 08

September 22, 2025

Review of Last Class

Dissolved Oxygen

Dissolved oxygen (DO) is the free, non-compound oxygen present in water or other liquids.

Freshwater can only hold small amounts, and this capacity is regulated by temperature.

Source: fondriest.com

Dissolved Oxygen and Life

Dissolved oxygen is an important nutrient for aquatic life.

Hypoxia occurs when DO levels are \(< 2\) mg/L.

Minimum DO requirements for freshwater fish

Source: fondriest.com

Impact of Paris On Seine DO, 1874

Dissolved Oxygen Downstream of Paris, 1874

Source: Dmitrieva, T., et al. (2018). https://doi.org/10.1007/s12685-018-0216-7

Modeling Dissolved Oxygen

Oxygen Balance in Rivers and Streams

Processes influencing oxygen balance in moving freshwater

DO Mass Balance

Let \(U\) be the river velocity (km/d), \(x\) the distance downstream from a waste release site in km, and \(C(x)\) the DO concentration at \(x\) in mg/L.

\[\begin{align} U \frac{dC}{dx} &= \Delta\ \text{DO} \\[0.5em] &= \text{Reaeration} + \text{Photosynthesis} - \text{Respiration} \\[0.5em] & \qquad - \text{Benthal Uptake} - \text{CBOD} - \text{NBOD} \end{align}\]

DO Mass Balance

Putting it all together:

\[ \begin{aligned} U \frac{dC}{dx} &= k_a (C_s - C) + P - R - S_B \\ &\quad - k_cB_0\exp\left(\frac{-k_cx}{U}\right) - k_n N_0\exp\left(\frac{-k_nx}{U}\right) \end{aligned} \]

Note: Usually models ignore \(P\), \(R\), and \(S_B\).

Why do you think that might be?

Analytic Solution (Streeter-Phelps)

\[\begin{align} C(x) &= C_s(1 - \alpha_1) + C_0 \alpha_1 - B_0 \alpha_2 - N_0 \alpha_3 + \left(\frac{P-R-S_B}{k_a}\right) (1-\alpha_1), \\[1em] \alpha_1 &= \exp\left(-\frac{k_a x}{U}\right) \\[0.25em] \alpha_2 &= \left(\frac{k_c}{k_a-k_c}\right)\left[\exp\left(\frac{-k_c x}{U}\right) - \exp\left(\frac{-k_ax}{U}\right)\right] \\[0.25em] \alpha_3 &= \left(\frac{k_n}{k_a-k_n}\right)\left[\exp\left(\frac{-k_n x}{U}\right) - \exp\left(\frac{-k_ax}{U}\right)\right] \end{align}\]

Simulation Outputs (“Sag Curve”)

Code

function do_analytic(x, C0, B0, N0, ka, kn, kc, Cs, U)
    B = B0 * exp(-kc * x / U)
    N = N0 * exp(-kn * x / U)
    α1 = exp(-ka * x / U)
    α2 = (kc/(ka-kc)) * (exp.(-kc * x / U) - exp(-ka * x / U))
    α3 = (kn/(ka-kn)) * (exp(-kn * x / U) - exp(-ka * x / U))
    C = Cs * (1 - α1) + (C0 * α1) - (B0 * α2) - (N0 * α3)
    return (C, B, N)
end  

# set river properties
ka = 0.6
kc = 0.4
kn = 0.25

C0 = 6.2
B0 = 9
N0 = 7

Cs = 7
U = 5

x = 0:40

# evaluate model over all x's
# this uses broadcasting
do_out = (y -> do_analytic(y, C0, B0, N0, ka, kc, kn, Cs, U)).(x)
# unpack outputs into individual arrays for C, B, and N
# this uses comprehensions to pull out the relevant components 
#of the tuples that our function outputs
C = [d[1] for d in do_out]
B = [d[2] for d in do_out]
N = [d[3] for d in do_out]

# plot outputs
p1 = plot(; ylabel="DO/OD (mg/l)", xlabel="Distance (km)", left_margin=8mm, legend=:outerright, bottom_margin=10mm)
plot!(p1, x, C, color=:black, linewidth=4, label="DO")
plot!(p1, x, B, color=:green, label="CBOD", linestyle=:dash, linewidth=3)
plot!(p1, x, N, color=:blue, label="NBOD", linestyle=:dash, linewidth=3)
# plot Cs, which is a constant value 
plot!(p1, x, Cs * ones(length(x)), color=:purple, label=L"C_s", linestyle=:dot, linewidth=2)

hline!([3], color=:red, linewidth=2, label="Regulatory Standard")
plot!(size=(1200, 450))
xaxis!((0, 40))

Figure 1: Sag curve for dissolved oxygen

Questions

Text: VSRIKRISH to 22333

URL: https://pollev.com/vsrikrish

See Results

Multiple Waste Discharges Impact on DO

Multiple Discharges

What happens if we have multiple discharge sites?

Schematic for Multiple Discharge Example

Simulating Multiple Discharges

Think of this as a multi-box system.

Flow from waste release 1 to waste release 2.
Flow from waste release 2 on.

How do we compute the initial conditions at release 2?

Who is Responsible for Non-Compliance?

Figure 2: Multi-Release Dissolved Oxygen Example

Numerical Simulation of DO ODE

Numerical Solution

What is the discretized version of the DO ODE?

\[ \begin{aligned} U \frac{dC}{dx} &= k_a (C_s - C) + P - R - S_B \\ &\quad - k_cB_0\exp\left(\frac{-k_cx}{U}\right) - k_n N_0\exp\left(\frac{-k_nx}{U}\right) \end{aligned} \]

Simulation Model Workflow

Numerical Solution

Divide river length into discrete spatial boxes of resolution \(\Delta x\) (e.g. 1km resolution).

\[ \begin{aligned} C(x+\Delta x) &= C(x) + \frac{\Delta x}{U}(k_a (C_s - C(x)) + P - R - S_B \\ &\quad - k_cB(x) - k_n N(x)) \\[0.5em] B(x + \Delta x) &= B(x)\exp\left(\frac{-k_c\Delta x}{U}\right) \\ N(x + \Delta x) &= N(x)\exp\left(\frac{-k_n\Delta x}{U}\right) \end{aligned} \]

Numerical Implementation

function do_numerical(L, Δx, C0, B0, N0, ka, kn, kc, Cs, U)
    n = Int(L / Δx) # number of length steps
    C = zeros(n + 1)
    B = zeros(n + 1)
    N = zeros(n + 1)
    C[1] = C0
    B[1] = B0
    N[1] = N0
    for i = 1:n
        B[i+1] = B[i] * exp(-kc * Δx / U)
        N[i+1] = N[i] * exp(-kn * Δx / U)
        C[i+1] = C[i] + (Δx / U) * (ka * (Cs - C[i]) - kc * B[i] - kn * N[i])
    end
    return (C, B, N)
end

Numerical vs. Analytic Solution

Code

C1, B1, N1 = do_numerical(40, 1, C0, B0, N0, ka, kc, kn, Cs, U)
C2, B2, N2 = do_numerical(40, 0.5, C0, B0, N0, ka, kc, kn, Cs, U)
C3, B3, N3 = do_numerical(40, 0.25, C0, B0, N0, ka, kc, kn, Cs, U)

p2 = plot(; ylabel="DO/OD (mg/l)", xlabel="Distance (km)", left_margin=10mm)
plot!(p2, 0:1:40, C, color=:black, linewidth=4, label="Analytic")
plot!(p2, 0:1:40, C1, color=:blue, linewidth=3, linestyle=:dash, label="Numerical (Δx = 1km)")
plot!(p2, 0:0.5:40, C2, color=:orange, linewidth=3, linestyle=:dash, label="Numerical (Δx = 0.5km)")
plot!(p2, 0:0.25:40, C3, color=:purple, linewidth=3, linestyle=:dash, label="Numerical (Δx = 0.25km)")
ylims!(p2, (2.25, 6.25))
plot!(size=(1200, 500))

Figure 3: Numerical vs. Analytic Solutions for DO model

Pros and Cons of Numerical Simulation

Defined more generally (analytic solution undefined when \(k_a = k_c\) or \(k_a = k_n\))
Can add more complexity (time/temperature-variation by modeling change in time)
But be careful of numerical errors and computational complexity.

Validation and Calibration of Simulation Models

Calibration vs. Validation

Framework from Oreskes et al. (1994):

Calibration: Selection of model parameters and structures to match simulations with observations.
Validation: Establishment of model legitimacy through evidence of consistency between simulation and real-world dynamics.

How Do We Calibrate?

Typical procedures include:

Selection of parameters based on experimental results;
“Curve-fitting”: minimize some error metric, such as root-mean-square-error: \[\text{RMSE} = \sqrt{\frac{1}{n}\sum_{i=1^n} (y_\text{sim} - y_\text{obs})^2}\]

Curve-fitting can run into problems when systems have multiple compensating processes.

Nonuniqueness: multiple parameter combinations/model structures provide equally good fits to observations.

Calibration Example

Suppose we don’t know \(k_a\), \(k_c\), and \(k_n\) but we have data:

Code

sim_dat = (x -> do_analytic(x, C0, B0, N0, 0.6, 0.4, 0.3, Cs, U)).(0:40)
C_dat = [d[1] for d in sim_dat]
C_noise = C_dat + rand(Normal(0, 0.4), length(C_dat))
p_calib = scatter(0:40, C_noise, xlabel="Distance (km)", ylabel="Dissolved Oxygen (mg/L)", label="Data", color=:blue, markersize=5)
plot!(p_calib, 0:40, C_dat, label="True Model", linewidth=3, color=:blue)
plot!(p_calib, size=(1000, 400))

Figure 4: Calibration data example with DO.

RMSE Minimization

rmse(x, y) = sqrt(mean((y .- x).^2))

function do_rmse(p, dat)
    ka, kc, kn = p # unpack parameter vector into unknown model variables
    C_sim, B, N = do_numerical(40, 0.25, C0, B0, N0, ka, kc, kn, Cs, U)
    return rmse(C_sim[1:4:length(C_sim)], dat)
end

lb = zeros(3)
ub = ones(3)
x0 = 0.5 * ones(3)
optim_out = Optim.optimize(p -> do_rmse(p, C_noise), lb, ub, x0)
@show round.(optim_out.minimizer; digits=1);
@show round.(optim_out.minimum; digits=2);

round.(optim_out.minimizer; digits = 1) = [0.6, 0.3, 0.3]
round.(optim_out.minimum; digits = 2) = 0.43

Calibration Results

Code

p_fit = optim_out.minimizer
C_fit, B, N = do_numerical(40, 0.25, C0, B0, N0, p_fit[1], p_fit[2], p_fit[3], Cs, U)
plot!(p_calib, 0:0.25:40, C_fit, linewidth=3, color=:orange, label="Calibrated Model")
plot!(p_calib, size=(1000, 500))

Figure 5: Calibration data example with DO.

Nonuniqueness

The calibration model produced the “wrong” parameter values:

Fitted: [0.6, 0.3, 0.3]
True: \(0.6, 0.4, 0.3\).

The RMSE was slightly lower than for the “true” model (0.427 vs 0.429)…

What can we attribute this “error” to? Do we care? How could we fix it?

Positive and Negative Controls

What we just did is called a “positive control” experiment: we tested the ability of our numerical model + calibration procedure to recover the parameters for synthetically generated data (so we know the “true” parameters).

We can also do a negative control: if we feed the model/procedure no signal (e.g. white noise), do we falsely recover a signal?

How Do We Validate?

Typical procedures include:

Compare simulations to out-of-sample data (not used for calibration!) (empirical validity).
Mechanistic justifications/assessments about which processes were included/not included and how they were parameterized (face validity).

Key Points

Attributing blame/liability for environmental degradation can be challenged by system complexity.
Numerical simulation example: be careful of potential for numerical errors, particularly near thresholds.
Calibration vs. validation: language often used in very loose fashion, but “calibration” refers to parameter selection and “validation” to assessment of model adequacy.

Upcoming Schedule and References

Next Classes

Wednesday: Uncertainty and Risk

Next Week: Monte Carlo Simulation

Assessments

Homework 2: Due Thursday (9/25) at 9pm.

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