Puff GaussianKC Model

This document demonstrates the setup and execution of a Puff Model for simulating the trajectory and concentration of multiple puffs. The model computes how puffs disperse in the atmosphere over time. This example uses GaussianKC model that calculates the time evolution of horizontal puff dispersion (sigmax, sigmay) based on turbulent velocity fluctuations (σux, σuy). It also computes the ground-level centerline concentration (Cgl) per unit mass assuming a top-hat vertical distribution within the lowest layer. The concentration is evaluated only when the puff is within the surface layer (zagl ≤ Δz); otherwise, it is set to zero. GaussianKC model must be coupled with BoundaryLayerMixingKC model that implements the Kantha-Clayson & Garratt turbulence parameterization for boundary layer mixing, as described in NOAA ARL-224 (https://www.arl.noaa.gov/documents/reports/arl-224.pdf) and HYSPLIT documentation.

We begin by importing the necessary libraries.

using Dates
using Plots
using Random
using EarthSciMLBase, EarthSciData, EnvironmentalTransport
using Aerosol
using ModelingToolkit, StochasticDiffEq
using ModelingToolkit: t
using EnvironmentalTransport: GaussianKC, PuffCoupler, GaussianKCCoupler, BoundaryLayerMixingKC, BoundaryLayerMixingKCCoupler

Reproducibility Setup

Here we set a fixed seed (e.g., 1234) so the random turbulence is reproducible across different runs.

Random.seed!(1234)

Random.TaskLocalRNG()

Simulation Parameters

The following are key simulation parameters used to define the puff model and the scenario being simulated:

start: Defines the time when the event starts. simulationlength: Duration of the simulation (24 hours in this case). puffreleaseinterval: The time interval between each puff release (in seconds). puffreleaseperinterval: The number of puffs released per time interval. ec0: Initial mass of elemental carbon released in each interval.

start = DateTime(2019, 06, 15, 0, 0, 0)
simulationlength = 24 * 1
puff_release_interval = 3600.0
puff_release_per_interval = 5
ec_0 = 100.0

puffs = [
    (-120.6, 46.6, 3.0),
    (-119.0, 47.0, 5.0),
    (-121.0, 46.5, 4.0),
    (-120.8, 46.8, 3.5),
    (-120.4, 46.7, 6.0),
    (-120.9, 46.6, 2.0),
    (-119.5, 46.9, 3.0),
    (-120.7, 47.2, 5.5),
    (-120.2, 46.5, 4.0),
    (-120.6, 46.3, 6.0)
]

Δλ_deg = 0.3125
Δφ_deg = 0.25
Δλ = deg2rad(Δλ_deg)
Δφ = deg2rad(Δφ_deg)

0.004363323129985824

Domain and Model

We define the simulation domain, including the range of longitudes, latitudes, and vertical levels. Then, we couple the models necessary for simulating particle dispersion and mixing, such as Puff, BoundaryLayerMixingKC, and GaussianKC.

DomainInfo: Defines the spatial and temporal extent of the simulation, including longitude, latitude, and vertical levels.

domain = DomainInfo(
    start, start + Hour(simulationlength);
    lonrange = deg2rad(-130 - Δλ_deg):Δλ:deg2rad(-60 + Δλ_deg),
    latrange = deg2rad(25 - Δφ_deg):Δφ:deg2rad(61 + Δφ_deg),
    levrange = 1:72
)

model = couple(
    Puff(domain),
    BoundaryLayerMixingKC(),
    GEOSFP("0.25x0.3125", domain; stream = false),
    ElementalCarbon(),
    GaussianKC()
)

const sys = convert(System, model)
equations(sys)[1:5]

\[ \begin{align} \frac{\mathrm{d} ~ \mathtt{Puff.lon}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{GEOSFP.A3dyn.U\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right), clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right)}{\mathtt{GEOSFP.lon2m} ~ \cos\left( \mathtt{Puff.lat}\left( t \right) \right)} + \frac{\mathtt{BoundaryLayerMixingKC.uprime\_x}\left( t \right)}{\mathtt{GEOSFP.lon2m} ~ \cos\left( \mathtt{Puff.lat}\left( t \right) \right)} \\ \frac{\mathrm{d} ~ \mathtt{Puff.lev}\left( t \right)}{\mathrm{d}t} &= ifelse\left( - \mathtt{Puff.offset} + \mathtt{Puff.lev}\left( t \right) < \mathtt{Puff.glo}, max\left( \mathtt{Puff.v\_zero}, \frac{\mathtt{GEOSFP.A3dyn.OMEGA\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right), clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right)}{\mathtt{GEOSFP.P\_unit} ~ derivative\left( Interpolation\_interp, clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right), 1 \right) + derivative\left( Interpolation\_interp, clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right), 1 \right) ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \right), ifelse\left( \mathtt{Puff.offset} + \mathtt{Puff.lev}\left( t \right) > \mathtt{Puff.ghi}, min\left( \mathtt{Puff.v\_zero}, \frac{\mathtt{GEOSFP.A3dyn.OMEGA\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right), clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right)}{\mathtt{GEOSFP.P\_unit} ~ derivative\left( Interpolation\_interp, clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right), 1 \right) + derivative\left( Interpolation\_interp, clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right), 1 \right) ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \right), \frac{\mathtt{GEOSFP.A3dyn.OMEGA\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right), clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right)}{\mathtt{GEOSFP.P\_unit} ~ derivative\left( Interpolation\_interp, clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right), 1 \right) + derivative\left( Interpolation\_interp, clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right), 1 \right) ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \right) \right) + \frac{ - \left( \mathtt{\_Pu\_blm} ~ DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64}([0.0, 4.804826, 659.3752000000001, 1313.48, 1961.311, 2609.201, 3257.081, 3898.201, 4533.901, 5169.611, 5805.321, 6436.264, 7062.197999999999, 7883.4220000000005, 8909.992, 9936.521, 10918.17, 11895.86, 12869.59, 14291.0, 15626.0, 16960.9, 18161.9, 19309.7, 20325.899999999998, 21215.0, 21877.600000000002, 22389.8, 22436.3, 21686.5, 20119.2, 17693.0, 15039.3, 12783.7, 10866.3, 9236.572, 7851.231, 6660.340999999999, 5638.791, 4764.391, 4017.541, 3381.0009999999997, 2836.781, 2373.0409999999997, 1979.1599999999999, 1645.71, 1364.34, 1127.69, 929.2942, 761.9842, 621.6801, 504.68010000000004, 407.6571, 327.6431, 262.0211, 208.497, 165.079, 130.05100000000002, 101.94399999999999, 79.51341, 61.677910000000004, 47.58061, 36.50411, 27.85261, 21.134900000000002, 15.9495, 11.9703, 8.934502, 6.600001, 4.758501, 3.27, 2.0, 1.0], 1:73, Float64[], DataInterpolations.LinearParameterCache{Vector{Float64}}(Float64[]), DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{UnitRange{Int64}}(1:73, Base.RefValue{Int64}(1), true), false, true)\left( 0.5 + clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right) + DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64}([1.0, 0.984952, 0.963406, 0.941865, 0.920387, 0.898908, 0.877429, 0.856018, 0.8346609, 0.8133039, 0.7919469, 0.7706375, 0.7493782, 0.721166, 0.6858999, 0.6506349, 0.6158184, 0.5810415, 0.5463042, 0.4945902, 0.4437402, 0.3928911, 0.3433811, 0.2944031, 0.2467411, 0.2003501, 0.1562241, 0.1136021, 0.06372006, 0.02801004, 0.006960025, 8.175413e-9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 1:73, Float64[], DataInterpolations.LinearParameterCache{Vector{Float64}}(Float64[]), DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{UnitRange{Int64}}(1:73, Base.RefValue{Int64}(1), true), false, true)\left( 0.5 + clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right) ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{\_g\_blm} ~ ifelse\left( \mathtt{Puff.lev}\left( t \right) \leq \mathtt{BoundaryLayerMixingKC.lev\_gnd} + \mathtt{BoundaryLayerMixingKC.offset}, \left( \mathtt{BoundaryLayerMixingKC.one} - \mathtt{BoundaryLayerMixingKC.w\_damp} \right) ~ \left|\mathtt{BoundaryLayerMixingKC.wprime}\left( t \right)\right|, ifelse\left( \mathtt{Puff.lev}\left( t \right) \geq \mathtt{BoundaryLayerMixingKC.lev\_top} - \mathtt{BoundaryLayerMixingKC.offset}, \left( - \mathtt{BoundaryLayerMixingKC.one} + \mathtt{BoundaryLayerMixingKC.w\_damp} \right) ~ \left|\mathtt{BoundaryLayerMixingKC.wprime}\left( t \right)\right|, \mathtt{BoundaryLayerMixingKC.wprime}\left( t \right) \right) \right)}{\left( - \mathtt{\_Pu\_blm} ~ DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64}([0.0, 4.804826, 659.3752000000001, 1313.48, 1961.311, 2609.201, 3257.081, 3898.201, 4533.901, 5169.611, 5805.321, 6436.264, 7062.197999999999, 7883.4220000000005, 8909.992, 9936.521, 10918.17, 11895.86, 12869.59, 14291.0, 15626.0, 16960.9, 18161.9, 19309.7, 20325.899999999998, 21215.0, 21877.600000000002, 22389.8, 22436.3, 21686.5, 20119.2, 17693.0, 15039.3, 12783.7, 10866.3, 9236.572, 7851.231, 6660.340999999999, 5638.791, 4764.391, 4017.541, 3381.0009999999997, 2836.781, 2373.0409999999997, 1979.1599999999999, 1645.71, 1364.34, 1127.69, 929.2942, 761.9842, 621.6801, 504.68010000000004, 407.6571, 327.6431, 262.0211, 208.497, 165.079, 130.05100000000002, 101.94399999999999, 79.51341, 61.677910000000004, 47.58061, 36.50411, 27.85261, 21.134900000000002, 15.9495, 11.9703, 8.934502, 6.600001, 4.758501, 3.27, 2.0, 1.0], 1:73, Float64[], DataInterpolations.LinearParameterCache{Vector{Float64}}(Float64[]), DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{UnitRange{Int64}}(1:73, Base.RefValue{Int64}(1), true), false, true)\left( clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right) + \mathtt{\_Pu\_blm} ~ DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64}([0.0, 4.804826, 659.3752000000001, 1313.48, 1961.311, 2609.201, 3257.081, 3898.201, 4533.901, 5169.611, 5805.321, 6436.264, 7062.197999999999, 7883.4220000000005, 8909.992, 9936.521, 10918.17, 11895.86, 12869.59, 14291.0, 15626.0, 16960.9, 18161.9, 19309.7, 20325.899999999998, 21215.0, 21877.600000000002, 22389.8, 22436.3, 21686.5, 20119.2, 17693.0, 15039.3, 12783.7, 10866.3, 9236.572, 7851.231, 6660.340999999999, 5638.791, 4764.391, 4017.541, 3381.0009999999997, 2836.781, 2373.0409999999997, 1979.1599999999999, 1645.71, 1364.34, 1127.69, 929.2942, 761.9842, 621.6801, 504.68010000000004, 407.6571, 327.6431, 262.0211, 208.497, 165.079, 130.05100000000002, 101.94399999999999, 79.51341, 61.677910000000004, 47.58061, 36.50411, 27.85261, 21.134900000000002, 15.9495, 11.9703, 8.934502, 6.600001, 4.758501, 3.27, 2.0, 1.0], 1:73, Float64[], DataInterpolations.LinearParameterCache{Vector{Float64}}(Float64[]), DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{UnitRange{Int64}}(1:73, Base.RefValue{Int64}(1), true), false, true)\left( 1 + clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right) + DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64}([1.0, 0.984952, 0.963406, 0.941865, 0.920387, 0.898908, 0.877429, 0.856018, 0.8346609, 0.8133039, 0.7919469, 0.7706375, 0.7493782, 0.721166, 0.6858999, 0.6506349, 0.6158184, 0.5810415, 0.5463042, 0.4945902, 0.4437402, 0.3928911, 0.3433811, 0.2944031, 0.2467411, 0.2003501, 0.1562241, 0.1136021, 0.06372006, 0.02801004, 0.006960025, 8.175413e-9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 1:73, Float64[], DataInterpolations.LinearParameterCache{Vector{Float64}}(Float64[]), DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{UnitRange{Int64}}(1:73, Base.RefValue{Int64}(1), true), false, true)\left( 1 + clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right) ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) - DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64}([1.0, 0.984952, 0.963406, 0.941865, 0.920387, 0.898908, 0.877429, 0.856018, 0.8346609, 0.8133039, 0.7919469, 0.7706375, 0.7493782, 0.721166, 0.6858999, 0.6506349, 0.6158184, 0.5810415, 0.5463042, 0.4945902, 0.4437402, 0.3928911, 0.3433811, 0.2944031, 0.2467411, 0.2003501, 0.1562241, 0.1136021, 0.06372006, 0.02801004, 0.006960025, 8.175413e-9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 1:73, Float64[], DataInterpolations.LinearParameterCache{Vector{Float64}}(Float64[]), DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{UnitRange{Int64}}(1:73, Base.RefValue{Int64}(1), true), false, true)\left( clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right) ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{\_Rd\_blm} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \\ \frac{\mathrm{d} ~ \mathtt{Puff.lat}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{BoundaryLayerMixingKC.uprime\_y}\left( t \right)}{\mathtt{GEOSFP.lat2meters}} + \frac{\mathtt{GEOSFP.A3dyn.V\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right), clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right)}{\mathtt{GEOSFP.lat2meters}} \\ \frac{\mathrm{d} ~ \mathtt{GaussianKC.sigma\_y}\left( t \right)}{\mathrm{d}t} &= 1.4142 ~ \sqrt{ifelse\left( \frac{ - \left( \mathtt{GEOSFP.A1.USTAR\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{3} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.g} ~ \mathtt{BoundaryLayerMixingKC.\kappa} ~ ifelse\left( \frac{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \mathtt{GEOSFP.A1.HFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.cp} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} + \frac{0.61 ~ \left( \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.EFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.Lv} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \geq \mathtt{BoundaryLayerMixingKC.kms0}, max\left( \frac{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \mathtt{GEOSFP.A1.HFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.cp} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} + \frac{0.61 ~ \left( \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.EFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.Lv} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}, \mathtt{BoundaryLayerMixingKC.kms1e\_8} \right), min\left( \frac{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \mathtt{GEOSFP.A1.HFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.cp} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} + \frac{0.61 ~ \left( \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.EFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.Lv} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}, - \mathtt{BoundaryLayerMixingKC.kms1e\_8} \right) \right)} < \mathtt{BoundaryLayerMixingKC.m0}, 0.36 ~ \left( ifelse\left( \frac{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \mathtt{GEOSFP.A1.HFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.cp} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} + \frac{0.61 ~ \left( \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.EFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.Lv} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} > \mathtt{BoundaryLayerMixingKC.kms0}, \left( \frac{\mathtt{BoundaryLayerMixingKC.g} ~ \left( \frac{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \mathtt{GEOSFP.A1.HFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.cp} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} + \frac{0.61 ~ \left( \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.EFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.Lv} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \right) ~ max\left( \mathtt{GEOSFP.A1.PBLH\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right), \mathtt{BoundaryLayerMixingKC.KMIX0} \right)}{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \right)^{0.33333}, \mathtt{BoundaryLayerMixingKC.ms0} \right) \right)^{2}, 5 ~ \left( \mathtt{GEOSFP.A1.USTAR\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \left( max\left( \mathtt{BoundaryLayerMixingKC.one} - max\left( min\left( \frac{\mathtt{\_Rd\_blm} ~ \log\left( \frac{\mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{\_Pu\_blm} ~ DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64}([0.0, 4.804826, 659.3752000000001, 1313.48, 1961.311, 2609.201, 3257.081, 3898.201, 4533.901, 5169.611, 5805.321, 6436.264, 7062.197999999999, 7883.4220000000005, 8909.992, 9936.521, 10918.17, 11895.86, 12869.59, 14291.0, 15626.0, 16960.9, 18161.9, 19309.7, 20325.899999999998, 21215.0, 21877.600000000002, 22389.8, 22436.3, 21686.5, 20119.2, 17693.0, 15039.3, 12783.7, 10866.3, 9236.572, 7851.231, 6660.340999999999, 5638.791, 4764.391, 4017.541, 3381.0009999999997, 2836.781, 2373.0409999999997, 1979.1599999999999, 1645.71, 1364.34, 1127.69, 929.2942, 761.9842, 621.6801, 504.68010000000004, 407.6571, 327.6431, 262.0211, 208.497, 165.079, 130.05100000000002, 101.94399999999999, 79.51341, 61.677910000000004, 47.58061, 36.50411, 27.85261, 21.134900000000002, 15.9495, 11.9703, 8.934502, 6.600001, 4.758501, 3.27, 2.0, 1.0], 1:73, Float64[], DataInterpolations.LinearParameterCache{Vector{Float64}}(Float64[]), DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{UnitRange{Int64}}(1:73, Base.RefValue{Int64}(1), true), false, true)\left( 0.5 + clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right) + DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64}([1.0, 0.984952, 0.963406, 0.941865, 0.920387, 0.898908, 0.877429, 0.856018, 0.8346609, 0.8133039, 0.7919469, 0.7706375, 0.7493782, 0.721166, 0.6858999, 0.6506349, 0.6158184, 0.5810415, 0.5463042, 0.4945902, 0.4437402, 0.3928911, 0.3433811, 0.2944031, 0.2467411, 0.2003501, 0.1562241, 0.1136021, 0.06372006, 0.02801004, 0.006960025, 8.175413e-9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 1:73, Float64[], DataInterpolations.LinearParameterCache{Vector{Float64}}(Float64[]), DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{UnitRange{Int64}}(1:73, Base.RefValue{Int64}(1), true), false, true)\left( 0.5 + clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right) ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \right) ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{\_g\_blm} ~ max\left( max\left( \mathtt{GEOSFP.A1.PBLH\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right), \mathtt{BoundaryLayerMixingKC.KMIX0} \right), \mathtt{BoundaryLayerMixingKC.m1} \right)}, 1 \right), 0 \right), 0 \right) \right)^{1.5} \right)} \\ \frac{\mathrm{d} ~ \mathtt{GaussianKC.sigma\_x}\left( t \right)}{\mathrm{d}t} &= 1.4142 ~ \sqrt{ifelse\left( \frac{ - \left( \mathtt{GEOSFP.A1.USTAR\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{3} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.g} ~ \mathtt{BoundaryLayerMixingKC.\kappa} ~ ifelse\left( \frac{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \mathtt{GEOSFP.A1.HFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.cp} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} + \frac{0.61 ~ \left( \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.EFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.Lv} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \geq \mathtt{BoundaryLayerMixingKC.kms0}, max\left( \frac{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \mathtt{GEOSFP.A1.HFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.cp} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} + \frac{0.61 ~ \left( \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.EFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.Lv} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}, \mathtt{BoundaryLayerMixingKC.kms1e\_8} \right), min\left( \frac{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \mathtt{GEOSFP.A1.HFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.cp} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} + \frac{0.61 ~ \left( \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.EFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.Lv} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}, - \mathtt{BoundaryLayerMixingKC.kms1e\_8} \right) \right)} < \mathtt{BoundaryLayerMixingKC.m0}, 0.36 ~ \left( ifelse\left( \frac{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \mathtt{GEOSFP.A1.HFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.cp} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} + \frac{0.61 ~ \left( \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.EFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.Lv} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} > \mathtt{BoundaryLayerMixingKC.kms0}, \left( \frac{\mathtt{BoundaryLayerMixingKC.g} ~ \left( \frac{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \mathtt{GEOSFP.A1.HFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.cp} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} + \frac{0.61 ~ \left( \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \mathtt{BoundaryLayerMixingKC.Rd} ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.EFLUX\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{BoundaryLayerMixingKC.Lv} ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \right) ~ max\left( \mathtt{GEOSFP.A1.PBLH\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right), \mathtt{BoundaryLayerMixingKC.KMIX0} \right)}{\left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \right)^{0.33333}, \mathtt{BoundaryLayerMixingKC.ms0} \right) \right)^{2}, 4 ~ \left( \mathtt{GEOSFP.A1.USTAR\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right)^{2} ~ \left( max\left( \mathtt{BoundaryLayerMixingKC.one} - max\left( min\left( \frac{\mathtt{\_Rd\_blm} ~ \log\left( \frac{\mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{\_Pu\_blm} ~ DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64}([0.0, 4.804826, 659.3752000000001, 1313.48, 1961.311, 2609.201, 3257.081, 3898.201, 4533.901, 5169.611, 5805.321, 6436.264, 7062.197999999999, 7883.4220000000005, 8909.992, 9936.521, 10918.17, 11895.86, 12869.59, 14291.0, 15626.0, 16960.9, 18161.9, 19309.7, 20325.899999999998, 21215.0, 21877.600000000002, 22389.8, 22436.3, 21686.5, 20119.2, 17693.0, 15039.3, 12783.7, 10866.3, 9236.572, 7851.231, 6660.340999999999, 5638.791, 4764.391, 4017.541, 3381.0009999999997, 2836.781, 2373.0409999999997, 1979.1599999999999, 1645.71, 1364.34, 1127.69, 929.2942, 761.9842, 621.6801, 504.68010000000004, 407.6571, 327.6431, 262.0211, 208.497, 165.079, 130.05100000000002, 101.94399999999999, 79.51341, 61.677910000000004, 47.58061, 36.50411, 27.85261, 21.134900000000002, 15.9495, 11.9703, 8.934502, 6.600001, 4.758501, 3.27, 2.0, 1.0], 1:73, Float64[], DataInterpolations.LinearParameterCache{Vector{Float64}}(Float64[]), DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{UnitRange{Int64}}(1:73, Base.RefValue{Int64}(1), true), false, true)\left( 0.5 + clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right) + DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Vector{Float64}, Vector{Float64}, Float64}([1.0, 0.984952, 0.963406, 0.941865, 0.920387, 0.898908, 0.877429, 0.856018, 0.8346609, 0.8133039, 0.7919469, 0.7706375, 0.7493782, 0.721166, 0.6858999, 0.6506349, 0.6158184, 0.5810415, 0.5463042, 0.4945902, 0.4437402, 0.3928911, 0.3433811, 0.2944031, 0.2467411, 0.2003501, 0.1562241, 0.1136021, 0.06372006, 0.02801004, 0.006960025, 8.175413e-9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 1:73, Float64[], DataInterpolations.LinearParameterCache{Vector{Float64}}(Float64[]), DataInterpolations.ExtrapolationType.None, DataInterpolations.ExtrapolationType.None, FindFirstFunctions.Guesser{UnitRange{Int64}}(1:73, Base.RefValue{Int64}(1), true), false, true)\left( 0.5 + clamp\left( \mathtt{Puff.lev}\left( t \right), 1, 72 \right) \right) ~ \mathtt{GEOSFP.I3.PS\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)} \right) ~ \left( 1 + 0.61 ~ \mathtt{GEOSFP.A1.QV2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right) \right) ~ \mathtt{GEOSFP.A1.T2M\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right)}{\mathtt{\_g\_blm} ~ max\left( max\left( \mathtt{GEOSFP.A1.PBLH\_itp}\left( \mathtt{GEOSFP.t\_ref} + t, \mathtt{Puff.lon}\left( t \right), \mathtt{Puff.lat}\left( t \right) \right), \mathtt{BoundaryLayerMixingKC.KMIX0} \right), \mathtt{BoundaryLayerMixingKC.m1} \right)}, 1 \right), 0 \right), 0 \right) \right)^{1.5} \right)} \end{align} \]

Solve The Model

In this step, we solve the model once to ensure that all necessary data is properly loaded and that the system is initialized correctly. We define the initial conditions for key parameters, such as sigmax (horizontal dispersion in the x-direction), sigmay (horizontal dispersion in the y-direction), and turbulent velocity components (wprime, uprimex, uprimey). We then set up the SDEProblem (Stochastic Differential Equation Problem) to model the system stochastically. The SRIW1 solver is used to solve the system of equations, accounting for the random variability in the system's behavior over time.

tspan = get_tspan(domain)

u0 = [
    sys.GaussianKC₊sigma_x => 0.00001,
    sys.GaussianKC₊sigma_y => 0.00001,
    sys.BoundaryLayerMixingKC₊wprime => 0.0,
    sys.BoundaryLayerMixingKC₊uprime_x => 0.0,
    sys.BoundaryLayerMixingKC₊uprime_y => 0.0,
]

prob = SDEProblem(sys, u0, tspan)

sol = solve(prob, SRIW1(); dt = 60.0)

retcode: Success
Interpolation: 1st order linear
t: 394-element Vector{Float64}:
     0.0
    60.0
   123.26084044251813
   194.42928594035104
   274.49378712541306
   364.56635095860787
   465.897985270952
   579.8960738723392
   708.1439235488998
   852.4227544350305
     ⋮
 80134.0437609088
 80830.15042477632
 81579.66198915709
 82422.86249908544
 83247.02221504935
 84174.20189550876
 85205.74746301718
 86320.87904533056
 86400.0
u: 394-element Vector{Vector{Float64}}:
 [-1.658062789394613, 36.5, 0.7504915783575618, 1.0e-5, 1.0e-5, 1.0, 0.0, 0.0, 0.0]
 [-1.6580461992773508, 36.49933330845963, 0.7504574080458638, 1.0e-5, 1.0e-5, 1.0, -0.016299325402272805, 0.0, 0.0]
 [-1.6580288665358036, 36.49804896102728, 0.7504214546781564, 1.0e-5, 1.0e-5, 1.0, -0.023001770296757106, 0.0, 0.0]
 [-1.6580095481386188, 36.49650619612452, 0.750381092019455, 1.0e-5, 1.0e-5, 1.0, -0.01847235460826546, 0.0, 0.0]
 [-1.6579880518239338, 36.494661150853496, 0.7503357950785295, 1.0e-5, 1.0e-5, 1.0, -0.021531597544144198, 0.0, 0.0]
 [-1.6579641530764773, 36.49290193146295, 0.7502849703420332, 1.0e-5, 1.0e-5, 1.0, -0.02353873703900995, 0.0, 0.0]
 [-1.6579376446308816, 36.49080006500533, 0.750227970544643, 1.0e-5, 1.0e-5, 1.0, -0.02614215923627773, 0.0, 0.0]
 [-1.6579082915576002, 36.488989400309094, 0.7501640679350594, 1.0e-5, 1.0e-5, 1.0, -0.008694023331678681, 0.0, 0.0]
 [-1.6578759549670194, 36.487105057574716, 0.7500924997977687, 1.0e-5, 1.0e-5, 1.0, -0.01719807677534944, 0.0, 0.0]
 [-1.657840367131711, 36.48379685025906, 0.7500123618704186, 1.0e-5, 1.0e-5, 1.0, -0.03059921544390249, 0.0, 0.0]
 ⋮
 [-1.645456857338743, 37.343698121753455, 0.7260356778028418, 22965.83969720238, 22965.83969720238, 1.0, 0.023844324754530726, -0.44558718657692914, -0.44558718657692914]
 [-1.6453245862401007, 37.36560977497005, 0.7258159558702999, 22965.83969720238, 22965.83969720238, 1.0, 0.003127752992058879, -0.4177927318128334, -0.4177927318128334]
 [-1.6451674179129656, 37.36148834270884, 0.7255716805495896, 22965.83969720238, 22965.83969720238, 1.0, -0.02607285562396005, -0.38980434000671205, -0.38980434000671205]
 [-1.6449808894273013, 37.34233715860447, 0.7252888568933807, 22965.83969720238, 22965.83969720238, 1.0, -0.010069244628157704, -0.36055874819730815, -0.36055874819730815]
 [-1.644789303379237, 37.33129397223276, 0.7250033415992457, 22965.83969720238, 22965.83969720238, 1.0, -0.02936941872856768, -0.33409395458582725, -0.33409395458582725]
 [-1.6445605451430627, 37.32395855996808, 0.7246709995386239, 22965.83969720238, 22965.83969720238, 1.0, -0.00021889144639912776, -0.30664317101148814, -0.30664317101148814]
 [-1.6442985131003742, 37.32408431173533, 0.724281453613618, 22965.83969720238, 22965.83969720238, 1.0, 0.02037450107362659, -0.27875334102555427, -0.27875334102555427]
 [-1.6440122710680112, 37.321802406718135, 0.7238347010155023, 22965.83969720238, 22965.83969720238, 1.0, -0.013249040423092526, -0.25145716128254086, -0.25145716128254086]
 [-1.643991511113184, 37.32117408313511, 0.7238020918728552, 22965.83969720238, 22965.83969720238, 1.0, 0.0018284076068899736, -0.24962173043181177, -0.24962173043181177]

Define Per-Puff Problems

We define the function prob_func, which sets the initial conditions for each puff. The initial release time for each puff is calculated based on the index i of the puff, and the release interval.

function prob_func(prob, i, repeat)
    rlon = deg2rad(puffs[i][1])
    rlat = deg2rad(puffs[i][2])
    rlev = puffs[i][3]

    u0 = [
        sys.Puff₊lon => rlon,
        sys.Puff₊lat => rlat,
        sys.Puff₊lev => rlev,
        sys.ElementalCarbon₊EC => ec_0 / puff_release_per_interval,
    ]

    p = [
        sys.GaussianKC₊Δz => 100,
    ]

    ts = (tspan[1] + floor((i-1) / puff_release_per_interval) * puff_release_interval, tspan[2])
    remake(prob; u0 = u0, p = p, tspan = ts)
end

prob_func (generic function with 1 method)

Solve Each Puff and Plot Trajectories

We solve the problem for puffs and store the results. Then, we generate an animation of the puff trajectories on a 3D map to visualize their movement and dispersion over time.

ensemble_prob = EnsembleProblem(prob, prob_func = prob_func, safetycopy = false)

esol = solve(ensemble_prob, SRIW1(), EnsembleThreads(); trajectories = length(puffs), dt = 60.0)

vars = [sys.Puff₊lon, sys.Puff₊lat, sys.Puff₊lev]
ranges = [(Inf, -Inf), (Inf, -Inf), (Inf, -Inf)]

for sol in esol
    for (i, var) in enumerate(vars)
        if !isempty(sol)
            data = sol[var]
            current_min, current_max = ranges[i]
            ranges[i] = (min(current_min, minimum(data)), max(current_max, maximum(data)))
        end
    end
end

t_ref = EarthSciMLBase.get_tref(domain)
sim_end_time = start + Hour(simulationlength)

anim = @animate for dt in datetime2unix(start):puff_release_interval:datetime2unix(sim_end_time)
    t = dt - t_ref

    p = plot(
        xlim = rad2deg.(ranges[1]),
        ylim = rad2deg.(ranges[2]),
        zlim = ranges[3],
        title = "Time: $(unix2datetime(dt))",
        xlabel = "Longitude", ylabel = "Latitude", zlabel = "Level",
        camera = (30, 45)
    )

    for sol in esol
        if t < sol.t[1] || t > sol.t[end]
            continue
        end

        lon = rad2deg(sol(t, idxs=sys.Puff₊lon))
        lat = rad2deg(sol(t, idxs=sys.Puff₊lat))
        lev = sol(t, idxs=sys.Puff₊lev)
        conc = sol(t, idxs=sys.GaussianKC₊C_gl)

        marker_color = cgrad(:inferno)[min(1.0, log10(conc + 1e-9) / 5)]

        scatter!(p,
            [lon], [lat], [lev],
            label = :none,
            markercolor = marker_color,
            markersize = 5
        )
    end
end

gif(anim, "puff_animation.gif", fps = 5)

Plot Ground-Level Centerline Concentration

Next we plot the Ground-Level Centerline Concentration for all simulated puffs. This represents the concentration of the pollutant at the surface (per unit mass) at the center of each puff derived from the Gaussian puff equations.

p_conc_0 = plot(
    title = "Ground-level centerline concentration",
    xlabel = "Time",
    ylabel = "Concentration (1/m³)",
    legend = :outertopright,
    xrotation = 45
)

for (i, sol) in enumerate(esol)
    time_vals = [start + Second(round(Int, t)) for t in sol.t]
    c_gl_vals = sol[sys.GaussianKC₊C_gl]

    plot!(p_conc_0, time_vals, c_gl_vals,
          label = "Puff $i", lw = 2, alpha = 0.8)
end

p_conc_0

Concentration at a Receptor Point

In this section, we compute the concentration at a receptor location (e.g., lon = -120°, lat = 46.5°) by summing the contribution from all simulated puffs at each time step. For each puff, we evaluate an elliptical Gaussian kernel using the puff’s horizontal spreads (σx, σy), its center position (lon, lat), and its ground-level centerline concentration (C_gl) per unit mass. Finally, we plot concentration versus time.

receptor_lon_deg = -120.0
receptor_lat_deg = 46.5
λr = deg2rad(receptor_lon_deg)
φr = deg2rad(receptor_lat_deg)

R_earth = 6.371e6 # Earth radius (m) for local tangent-plane distance approximation

# dx ~ east-west distance, dy ~ north-south distance
dxdy_m(lon_puff, lat_puff, lon_rec, lat_rec) = begin
    dx = R_earth * cos(0.5*(lat_puff + lat_rec)) * (lon_rec - lon_puff)
    dy = R_earth * (lat_rec - lat_puff)
    dx, dy
end

# Build a time grid for plotting (e.g., every 1 hour)
sim_end_time = start + Hour(simulationlength)
tgrid_unix = collect(datetime2unix(start):3600.0:datetime2unix(sim_end_time))
tgrid = tgrid_unix .- t_ref
time_dt = unix2datetime.(tgrid_unix)

C_rec = zeros(length(tgrid)) # Total concentration at receptor

for (k, tsec) in enumerate(tgrid)
    csum = 0.0
    for sol in esol
        if isempty(sol) || tsec < sol.t[1] || tsec > sol.t[end]
            continue
        end

        # Puff center location
        lonp = sol(tsec, idxs = sys.Puff₊lon)
        latp = sol(tsec, idxs = sys.Puff₊lat)

        # GaussianKC horizontal dispersion std dev
        sx = sol(tsec, idxs = sys.GaussianKC₊sigma_x)
        sy = sol(tsec, idxs = sys.GaussianKC₊sigma_y)

        # Centerline ground-level concentration per unit mass (1/m^3)
        Cgl = sol(tsec, idxs = sys.GaussianKC₊C_gl)

        # Mass of elemental carbon in the Puffs
        m = sol(tsec, idxs = sys.ElementalCarbon₊EC)

        dx, dy = dxdy_m(lonp, latp, λr, φr)
        kernel = exp(-0.5 * ((dx / sx)^2 + (dy / sy)^2))

        csum += m * Cgl * kernel
    end
    C_rec[k] = csum
end

# Convert to µg/m^3
C_rec_ugm3 = C_rec .* 1e9

p_conc = plot(
    time_dt, C_rec_ugm3;
    title  = "Concentration at receptor (lon=$(receptor_lon_deg), lat=$(receptor_lat_deg))",
    xlabel = "Time",
    ylabel = "Concentration (µg/m³)",
    lw = 2,
    legend = false,
    xrotation = 45
)

p_conc

Animated Spatial Heatmap

Here we compute concentration on a 2D longitude–latitude grid at each time step using the same Gaussian kernel method as the single-receptor calculation above, then animate the result as a heatmap.

# Define spatial grid
lon_range = range(-122.0, -118.0; length=50)
lat_range = range(45.5, 48.0; length=50)

anim_heatmap = @animate for (k, tsec) in enumerate(tgrid)
    C_grid = zeros(length(lat_range), length(lon_range))

    for (j, lon_deg) in enumerate(lon_range)
        for (i, lat_deg) in enumerate(lat_range)
            λ_rec = deg2rad(lon_deg)
            φ_rec = deg2rad(lat_deg)
            csum = 0.0
            for sol in esol
                if isempty(sol) || tsec < sol.t[1] || tsec > sol.t[end]
                    continue
                end
                lonp = sol(tsec, idxs = sys.Puff₊lon)
                latp = sol(tsec, idxs = sys.Puff₊lat)
                sx   = sol(tsec, idxs = sys.GaussianKC₊sigma_x)
                sy   = sol(tsec, idxs = sys.GaussianKC₊sigma_y)
                Cgl  = sol(tsec, idxs = sys.GaussianKC₊C_gl)
                m    = sol(tsec, idxs = sys.ElementalCarbon₊EC)

                dx, dy = dxdy_m(lonp, latp, λ_rec, φ_rec)
                kernel = exp(-0.5 * ((dx / sx)^2 + (dy / sy)^2))
                csum += m * Cgl * kernel
            end
            C_grid[i, j] = csum * 1e9  # µg/m³
        end
    end

    heatmap(
        collect(lon_range), collect(lat_range), C_grid;
        title  = "Concentration at $(time_dt[k])",
        xlabel = "Longitude (°)",
        ylabel = "Latitude (°)",
        colorbar_title = "Concentration (µg/m³)",
        clims  = (0, maximum(C_rec_ugm3)),
        color  = :inferno
    )
end

gif(anim_heatmap, "puff_heatmap.gif", fps = 3)