Simulation and Box Models

Lecture 06

September 16, 2024

Review of Last Class

Shallow Lake Model

Model of phosphorous cycling
Features complex dynamics: multiple equilibria, feedbacks, bifurcations

Feedbacks

Unstable equilibria can result from reinforcing (positive) feedback loops, where a shock to the system state gets amplified.

Feedback loops can also be dampening (negative), where a shock is weakened (stable equilibria).

Tipping Points and Bifurcations

Uncertainty about system dynamics can dramatically change equilibria locations and behavior;
“Shocks” (in this case, sedimentation/recycling disturbances or massive non-point source inflows) can irreversibly alter system outcomes.

Questions?

Text: VSRIKRISH to 22333

URL: https://pollev.com/vsrikrish

See Results

Simulating Systems

What is Simulation?

Simulation: evaluating a model to understand how a system might evolve under a particular set of conditions.

Think of simulation as data generation (or generative modeling).
The model represents a particular data-generating process.

Why Simulate Systems?

System involves complex, nonlinear dynamics that may not be analytically tractable.
Setting up and running a real-world experiment is not possible.
State depends on prior states or states of nearby locations, so need to iterate over multiple spatial or temporal steps.
Need to understand range of system performance across rarely-seen conditions.

Types of Simulation Models

Deterministic vs. Stochastic
Discrete vs. Continuous

Simulation Model Workflow

Simulation Workflow

Simulation Model Applications

Water balance/hydrologic flow models;
Climate models (ocean heat/CO\(_2\) uptake through “box” layers)
Airsheds
Epidemiology
Social science (agent-based models)

Example: Airsheds

What Is A Box Model?

Box models are a common building block of simulation models.

Box models are all about mass-balance (mass \(m\)), assume well-mixed within box.

Can be steady-state \((\dot{m} = 0)\) or not.

Airshed Model

Let’s look at a simple steady-state model of an airshed.

Variable	Meaning	Units
\(m\)	mass of some air pollutant	g
\(C\)	concentration in box	g/m\(^3\)
\(S, D\)	source, deposition rate within the box	g/s
\(u\)	wind speed	m/s
\(L, W, H\)	box dimensions	m

Selecting Box Dimensions

What is relevant for the box dimensions \(L\), \(W\), and \(H\)?

Primarily the assumption(s) about mixing:

Mixing height: is there an atmospheric inversion which limits mixing height?
Homogeneity of input/output flows and emissions.

Steady-State Airshed Box Model

Steady-state box ⇒ \(\dot{m} = 0\).

\[\begin{align} 0 &= m_\text{in} - m_\text{out} + S - D \\[0.5em] &\class{\fragment}{{} = (u WH) C_\text{in} - (u WH) C + S - D } \\[0.5em] \end{align}\]

Solving for \(C\): \[C = C_{in} + \frac{S-D}{uWH}\]

Decay Processes

Now let’s assume some process affecting \(m\) depends on time.

For example: let’s say we care about an air pollutant which has a first-order decay rate \(k\), so \(D(t) = D_0 - km(t)\).

\[ \Rightarrow \frac{dm}{dt} = m_\text{in} - m_\text{out} + S - D_0 - km \]

\[ \dot{m} = \frac{d(CV)}{dt} = \overbrace{(u WH) C_\text{in}}^{\text{inflow}} - \overbrace{(u WH) C}^{\text{outflow}} + \overbrace{S - D_0}^{\text{net emissions}} - \overbrace{kCV}^{\text{mass decay}} \]

Discretizing Continuous Models

We could analytically solve this particular differential equation, but in general this may not be possible without strong assumptions.

Instead, we can discretize these models using methods from CEE 3200!

Euler Discretization

Recall that \[\frac{df}{dt} = \lim_{\Delta t \to 0} \frac{\Delta f(t)}{\Delta t}.\]

So if we pick a sufficiently small step size \(\Delta t\), can use this as an approximation:

\[\frac{df}{dt} \approx \frac{\Delta f(t)}{\Delta t}\]

Discretizing the Airshed Model

\[ \begin{align*} \frac{d(CV)}{dt} &= (u WH) C_\text{in} - (u WH) C(t) + S - D_0 - kC(t)V \\[0.5em] \frac{d(C)}{dt} &= \frac{u}{L} C_\text{in} - \frac{u}{L} C(t) + \frac{S - D_0}{V} - kC(t) \\ \end{align*} \]

\[ \frac{C(t+1) - C(t)}{\Delta t} = \frac{u}{L} C_\text{in} - \frac{u}{L} C(t) + \frac{S - D_0}{V} - kC(t) \]

\[ \bbox[5px, border: red 5px solid]{C(t+1) = C(t) + \Delta t \left(\frac{u}{L} \left(C_\text{in} - C(t)\right) + \frac{S - D_0}{V} - kC(t)\right)} \]

How To Simulate?

Pick \(\Delta t\);
Starting at \(t=1\), iterate equation until end-time \(T\).

Simple Simulation Code

# this function computes the increment by which the concentration
# is updated at each step
function box_simulate_timestep(C, Ci, u, W, H, L, S, D, k)
    return ((u / L) * (Ci - C) + (S - D) / (W * H * L) - (k * C))
end

# this function loops over the timesteps to simulate
# the concentration series
function airshed_simulate(C₀, Ci, u, W, H, L, S, D, k, T, Δt)
    # initialize C storage
    # for code simplicity we make the array length T+1 so
    # index 1 is C₀
    steps = Int64(T / Δt)
    C = zeros(steps + 1)
    C[1] = C₀
    for t = 1:steps
        C[t+1] = C[t] + 
            Δt * box_simulate_timestep(C[t], Ci, u, W, H, L, S, D, k)
    end
    # the first element of C is the initial condition
    return C
end
Δt = 1.0
T = 10
C₀ = 0.1
k = 0.3
u = 2
L = 4
W = 4
H = 4
Ci = 0.2
S = 10
D = 13
C = airshed_simulate(C₀, Ci, u, W, H, L, S, D, k, T, Δt)

Simulation Results

Code

p = plot(; xlabel="Time", ylabel=L"Pollutant concentration (g/m$^3$)")

# find exact solution and plot for comparison
# need these substitutions for the exact solution; not critical otherwise
V = 4^3
P = (u / L) * Ci + (S - D) / V
l = (u / L) + k
# use steps of 0.01 to smooth the plotting
C₁ = C₀ .* exp.(-l * (0:0.01:T))
C₂ = (P / l) * (1 .- exp.(-l * (0:0.01:T)))
C_exact = C₁ .+ C₂
plot!(p, 0:0.01:T, C_exact, linewidth=3, color=:black, label="Exact Solution")

# add simulated solution
plot!(p, 0:Δt:T, C, linewidth=3, color=:blue, label="Simulated Solution (Δt = 1)")

Figure 1: Simulation results for airshed model.

Impact of Time Step Size

Code

Csmall = airshed_simulate(C₀, Ci, u, W, H, L, S, D, k, T, 0.1)
plot!(p, 0:0.1:T, Csmall, linewidth=3, color=:purple, label="Simulated Solution (Δt = 0.1)")
Csmaller = airshed_simulate(C₀, Ci, u, W, H, L, S, D, k, T, 0.01)
plot!(p, 0:0.01:T, Csmaller, linewidth=3, color=:orange, label="Simulated Solution (Δt = 0.01)")

Figure 2: Simulation results for airshed model with varying timesteps.

Time-Varying (Dynamic) Simulation

In our prior example, inflow conditions were static (often an assumption needed for analytic solutions).

The simulation framework lets us make these time-varying:

\[ \begin{aligned} C(t+1) &= C(t) + \\ & \qquad \Delta t \left(\frac{{\color{red}u(t)}}{L} \left({\color{red}C_\text{in}(t)} - C(t)\right) + {\color{red}S(t)} - {\color{red}D(t)} - kC(t)\right) \end{aligned} \]

Dynamic Simulation Code

function airshed_simulate_dynamic(C₀, Ci, u, W, H, L, S, D, k, T, Δt)
    # initialize C storage
    # for code simplicity we make the array length T+1 so
    # index 1 is C₀
    steps = Int64(T / Δt)
    C = zeros(steps + 1)
    C[1] = C₀
    for t = 1:steps
        C[t+1] = C[t] + Δt * box_simulate_timestep(C[t], Ci[t], u[t], W, H, L, S, D, k)
    end
    # the first element of C is the initial condition
    return C
end

Dynamic Inputs

Code

Δt = 0.1
steps = Int(T / Δt) + 1
u = rand(LogNormal(log(2), 0.05), steps)
Ci = rand(LogNormal(log(0.2), 0.1), steps)

pwind = plot(0:Δt:T, u; title="Wind Speed", size=(550, 450), xlabel="Time", ylabel="m/s")
pincoming = plot(0:Δt:T, Ci; title="Incoming Concentration", size=(550, 450), xlabel="Time", ylabel=L"$\mathrm{g/m}^3$")

display(pwind)
display(pincoming)

(a) Dynamic simulation results for airshed model.

(b)

Figure 3

Dynamic Simulation Results

Code

Cdynamic = airshed_simulate_dynamic(C₀, Ci, u, W, H, L, S, D, k, T, 0.1)
p = plot(; xlabel="Time", ylabel=L"Pollutant Concentration ($\mathrm{g/m}^3$)")
plot!(p, 0:0.1:T, Cdynamic, linewidth=3, color=:purple, label="Dynamic Simulated Solution")
plot!(p, 0:0.1:T, Csmall, linewidth=3, color=:red, label="Static Simulated Solution")

Figure 4: Dynamic simulation results for airshed model.

Multi-Box Simulation

More Complex Domains

We can use the single box simulation as a building block for more complex domains, possibly with different dynamics.

Two box airshed model

Multi-Box Simulation Approach

Tip: Use smaller functions as a building block!

Multi-Box Simulation Approach

# here we use our timestep function on each box individually,
# but update them "together"
function airshed_twobox_simulate(C1₀, C2₀, Ci, u, W1, W2, H1, H2, L1, L2, S1, D1, S2, D2, k, T, Δt)
    # initialize C storage
    # for code simplicity we make the array length T+1 so
    # index 1 is C₀
    steps = Int64(T / Δt)
    C1 = zeros(steps + 1)
    C2 = zeros(steps + 1)
    C1[1] = C1₀
    C2[1] = C2₀
    # for each time step, first we update box 1, then box 2
    for t = 1:steps
        C1[t+1] = C1[t] + Δt * box_simulate_timestep(C1[t], Ci, u, W1, H1, L1, S1, D1, k)
        C2[t+1] = C2[t] + Δt * box_simulate_timestep(C2[t], C1[t+1], u, W2, H2, L2, S2, D2, k)
    end
    return (C1, C2)
end

Two-Box Simulation Results

Code

Δt = 0.1
T = 10
C1₀ = 0.05
C2₀ = 0.1
k = 0.3
u = 2
W1 = 6
H1 = 6
L1 = 6
W2 = 4
H2 = 4
L2 = 4
Ci = 0.01
S1 = 30
D1 = 20
S2 = 10
D2 = 13
C = airshed_twobox_simulate(C1₀, C2₀, Ci, u, W1, W2, H1, H2, L1, L2, S1, D1, S2, D2, k, T, Δt)

p1 = plot(0:0.1:T, C[1], linewidth=3, color=:purple, xlabel="Time", ylabel=L"Concentration ($\mathrm{g/m}^3$)", title="Box 1", size=(600, 500))
p2 = plot(0:0.1:T, C[2], linewidth=3, color=:red, xlabel="Time", ylabel=L"Concentration ($\mathrm{g/m}^3$)", title="Box 2", size=(600, 500))
display(p1)
display(p2)

(a) Simulation results for two-box airshed model.

(b)

Figure 5

Key Takeaways

Simulation modeling involves using a model to generate “data” under certain conditions.
Simulations are the main approach to descriptive modeling.
Divide domain into spatial domains and/or temporal steps and iterate.
Discretize equations if needed.

Simulation and Box Models

Review of Last Class

Shallow Lake Model

Feedbacks

Tipping Points and Bifurcations

Questions?

Simulating Systems

What is Simulation?

Why Simulate Systems?

Types of Simulation Models

Simulation Model Workflow

Simulation Model Applications

Example: Airsheds

What Is A Box Model?

Airshed Model

Selecting Box Dimensions

Steady-State Airshed Box Model

Decay Processes

Discretizing Continuous Models

Euler Discretization

Discretizing the Airshed Model

How To Simulate?

Simple Simulation Code

Simulation Results

Impact of Time Step Size

Time-Varying (Dynamic) Simulation

Dynamic Simulation Code

Dynamic Inputs

Dynamic Simulation Results

Multi-Box Simulation

More Complex Domains

Multi-Box Simulation Approach

Multi-Box Simulation Approach

Two-Box Simulation Results

Key Takeaways

Key Takeaways

Upcoming Schedule

Next Classes

Assessments

References

References