Create an EarthSciMLBase.Operator that performs advection. Advection is performed using the given stencil operator (e.g. l94_stencil or ppm_stencil). p is an optional parameter set to be used by the stencil operator. bc_type is the boundary condition type, e.g. ZeroGradBC().

Wind field data will be added in automatically if available. Currently the only valid source of wind data is EarthSciData.GEOSFP.

An array with external indexing implemented for boundary conditions.

Constant boundary conditions.

An array with zero constant boundary conditions.

PBLMixingCallback <: EarthSciMLBase.Operator

A callback that applies planetary boundary layer (PBL) mixing to tracer fields at periodic intervals. PBL mixing is a discrete process that redistributes tracers vertically within each grid column.

Species-specific constant boundary conditions. Takes a dictionary mapping species names/indices to boundary values and a default value.

Examples:

• SpeciesConstantBC(Dict("O3" => 40.0), 0.0) sets O3 to 40.0 and others to 0.0
• SpeciesConstantBC(Dict(1 => 40.0), 0.0) sets species 1 to 40.0 and others to 0.0
• SpeciesConstantBC(Dict("O3" => 40.0, "NO2" => 10.0), 0.0) sets multiple species

Note: When using species names, they will be resolved to indices when the boundary condition is applied to a system with known species variables.

An array with species-specific constant boundary conditions.

Zero gradient boundary conditions.

An array with zero gradient boundary conditions.

EarthSciMLBase.init_callback(cb::PBLMixingCallback, csys::CoupledSystem, sys_mtk, coord_args, domain::DomainInfo, alg)

Initialize the PBL mixing callback.

BoundaryLayerMixingKC()

A ModelingToolkit System implementing 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.

GaussianKC()

Returns a ModelingToolkit.System 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.

Note: Must be coupled with BoundaryLayerMixingKC.

Example:

using Dates, EarthSciMLBase, EarthSciData, EnvironmentalTransport
using ModelingToolkit, StochasticDiffEq

t0 = DateTime(2022, 5, 1)
t1 = DateTime(2022, 5, 2)
Δλ      = deg2rad(5.0)
Δφ      = deg2rad(4.0)

dom = DomainInfo(t0, t1;
                    lonrange = deg2rad(-130):Δλ:deg2rad(-60),
                    latrange = deg2rad(25):Δφ:deg2rad(61)),
                    levrange=1:72

mdl = couple(Puff(dom),
             BoundaryLayerMixingKC(),
             GEOSFP("4x5", dom; stream=false),
             GaussianKC())

sys  = convert(System, mdl)
tspan = get_tspan(dom)

u0 = [sys.Puff₊lon => deg2rad(-105),
      sys.Puff₊lat => deg2rad(  38),
      sys.Puff₊lev => 2,
      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]

p = [sys.GaussianKC₊Δz => 100.0]

prob = SDEProblem(sys, u0, tspan, p)
sol = solve(prob, SRIW1(); dt=60.0)

sigma_x = sol[sys.GaussianKC₊sigma_x]
sigma_y = sol[sys.GaussianKC₊sigma_y]
C_gl    = sol[sys.GaussianKC₊C_gl]

GaussianPGB()

Return a ModelingToolkit.System implementing a classic Gaussian plume dispersion model, parameterized with Pasquill-Gifford-Briggs dispersion coefficients, following the formulations described in EPA guidance 402-R-00-004 §12.1.6 (https://19january2017snapshot.epa.gov/sites/production/files/2015-05/documents/402-r-00-004.pdf) and the MMGRMA document, Table 6-7 (https://www.epa.gov/sites/default/files/2020-10/documents/mmgrma_0.pdf). Good for: quasi-steady (piecewise-steady) applications where emissions and meteorology can be treated steady over each model time step; near-field ranges (typically ≤ 50 km); and cases where the effective plume height remains within the planetary boundary layer.

What this does:

• Stability classification: Uses near-surface meteorology (10 m wind speed, downward short-wave radiation, total cloud fraction, and the surface temperature lapse) to determine the Pasquill stability class.

• Dispersion coefficients: Maps the stability class to Briggs coefficients (Ay, Az, By, Bz) and evaluates the analytic expressions for the horizontal (sigmah ≡ sigmay) and vertical (sigma_z) dispersion parameters as functions of the down-wind distance 'x'.

• Hypsometric height: Computes puff height above ground (z_agl) from pressure, temperature, and humidity using the hypsometric equation with a virtual temperature layer mean.

• Ground-level concentration: Computes the Gaussian ground-level concentration at the puff center for one unit of puff mass: Cgl = 1 / ((2π)^{3/2} · σh² · σz) * exp(-zagl² / (2 σ_z²)).

Example:

using Dates, EarthSciMLBase, EarthSciData, EnvironmentalTransport
using ModelingToolkit, OrdinaryDiffEq

t0 = DateTime(2022, 5, 1)
t1 = DateTime(2022, 5, 2)
Δλ      = deg2rad(5.0)
Δφ      = deg2rad(4.0)

dom = DomainInfo(t0, t1; levrange=1:72,
                    lonrange = deg2rad(-130):Δλ:deg2rad(-60),
                    latrange = deg2rad(25):Δφ:deg2rad(61))

mdl = couple(Puff(dom),
             GEOSFP("4x5", dom; stream=false),
             GaussianPGB())

sys  = convert(System, mdl)
tspan = get_tspan(dom)

u0 = [sys.Puff₊lon => deg2rad(-105),
      sys.Puff₊lat => deg2rad(  38),
      sys.Puff₊lev => 2]

p  = [sys.GaussianPGB₊lon0 => deg2rad(-105),
      sys.GaussianPGB₊lat0 => deg2rad(  38)]

prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, Tsit5();)

sigma_h = sol[sys.GaussianPGB₊sigma_h]
sigma_z = sol[sys.GaussianPGB₊sigma_z]
C_gl    = sol[sys.GaussianPGB₊C_gl]

HeavisideBoundaryCondition(
;
    name
) -> ModelingToolkitBase.System

Create a smoothed Heaviside boundary condition component for coupling surface and subsurface water flow, following Wang et al. (2020).

The Heaviside step function determines whether the soil surface is in a ponding state (h ≥ 0) or unsaturated state (h < 0), enabling automatic switching between infiltration and evaporation boundary conditions:

\[\eta_\omega(h) = \frac{1}{\pi} \arctan\left(\frac{h}{\omega}\right) + \frac{1}{2}\]

\[\delta_\omega(h) = \frac{1}{\pi} \frac{\omega}{\omega^2 + h^2}\]

These smooth approximations converge to the Heaviside step function and Dirac delta function respectively as $\omega \to 0^+$.

The boundary condition equation is (Eq. 3):

\[\frac{\partial(\eta(h) \cdot h)}{\partial t} - (P - I) = 0\]

Reference: Wang, Z., Timlin, D., Kouznetsov, M., Fleisher, D., Li, S., Tully, K., & Reddy, V. (2020). Coupled model of surface runoff and surface-subsurface water movement. Advances in Water Resources, 137, 103499. https://doi.org/10.1016/j.advwatres.2019.103499

Puff(
    di::EarthSciMLBase.DomainInfo;
    buffer_cells,
    name
) -> ModelingToolkitBase.System

Create a Lagrangian transport model which advects a "puff" or particle of matter within a fluid velocity field.

Model boundaries are set by the DomainInfo argument. The model sets boundaries at the ground and model bottom and top, preventing the puff from crossing those boundaries. If the puff reaches one of the horizontal boundaries, the simulation is stopped.

Keyword arguments

• buffer_cells: The distance (expressed in a number of DomainInfo grid cells) to use as a buffer around the horizontal edge of the domain to avoid data loader interpolation errors. The effective size of the domain will be reduce by 2× this amount (default = 1)

SaintVenantPDE(
    L_domain,
    T_end;
    P_val,
    I_val,
    S_0_val,
    n_manning_val,
    g_val,
    h_min_val,
    h_init_val,
    q_init_val,
    name
) -> ModelingToolkitBase.PDESystem

Create a PDESystem for the Saint-Venant equations (Eq. 1 from Wang et al., 2020) suitable for spatial discretization with MethodOfLines.jl.

The system implements the full 1D Saint-Venant equations with proper SI units:

• Mass conservation (Eq. 1a): $\partial \tilde{h}/\partial t = -\partial q/\partial l + (P - I)$
• Momentum conservation (Eq. 1b): $\partial q/\partial t = -\partial F/\partial l + g \tilde{h}(S_0 - S_f)$
• Momentum flux (auxiliary): $F = q^2/\tilde{h} + g \tilde{h}^2/2$

where $S_f$ is Manning's friction slope computed with non-dimensionalized fractional exponents: $S_f = ((n/n_{ref})(q/q_{ref}))^2 / (\tilde{h}/h_{ref})^{10/3}$.

Arguments

• L_domain: Length of the spatial domain (m)
• T_end: Duration of the simulation (s)

Keyword Arguments

• P_val: Precipitation rate (m/s), default 70 mm/hr
• I_val: Infiltration rate (m/s), default 0.0
• S_0_val: Surface slope (dimensionless), default 0.01
• n_manning_val: Manning roughness coefficient (m^(-1/3)·s), default 0.03
• g_val: Gravitational acceleration (m/s²), default 9.81
• h_min_val: Minimum flow depth to prevent singularity (m), default 1e-5
• h_init_val: Initial and boundary flow depth (m), default 1e-3
• q_init_val: Initial and boundary flux (m²/s), default 0.0
• name: System name, default :SaintVenantPDE

Reference: Wang, Z., Timlin, D., Kouznetsov, M., Fleisher, D., Li, S., Tully, K., & Reddy, V. (2020). Coupled model of surface runoff and surface-subsurface water movement. Advances in Water Resources, 137, 103499. https://doi.org/10.1016/j.advwatres.2019.103499

Wildfire plume rise model based on Sofiev et al. (2012) [1].

[1] Sofiev, M., Ermakova, T., and Vankevich, R.: Evaluation of the smoke-injection height from wild-land fires using remote-sensing data, Atmos. Chem. Phys., 12, 1995–2006, https://doi.org/10.5194/acp-12-1995-2012, 2012.

SurfaceRunoff(; name) -> ModelingToolkitBase.System

Create a surface runoff model based on the Saint-Venant equation system for surface water movement, following the approach of Wang et al. (2020).

The model describes mass conservation and momentum conservation of surface water flow along a soil surface slope. The Saint-Venant equations are:

\[\frac{\partial \tilde{h}}{\partial t} = -\frac{\partial q}{\partial l} + (P - I)\]

\[\frac{\partial q}{\partial t} = -\frac{\partial}{\partial l}\left(\frac{q^2}{\tilde{h}} + \frac{g \tilde{h}^2}{2}\right) + g \tilde{h} (S_0 - S_f)\]

where $\tilde{h}$ is the flow depth (ponded water height), $q$ is the surface runoff flux per unit width, $P$ is precipitation flux, $I$ is infiltration flux, $S_0$ is the surface slope, and $S_f$ is the friction slope computed from Manning's equation.

This component represents the equations at a single surface node. The spatial derivative terms ($\partial q / \partial l$ and the momentum flux divergence) are represented as input parameters that should be provided by a spatial discretization scheme or coupled model.

Reference: Wang, Z., Timlin, D., Kouznetsov, M., Fleisher, D., Li, S., Tully, K., & Reddy, V. (2020). Coupled model of surface runoff and surface-subsurface water movement. Advances in Water Resources, 137, 103499. https://doi.org/10.1016/j.advwatres.2019.103499

An advection kernel for a 4D array, where the first dimension is the state variables and the next three dimensions are the spatial dimensions.

A function to create an advection operator for a 4D array,

Arguments:

• u_prototype: A prototype array of the same size and type as the input array.
• stencil: The stencil operator, e.g. l94_stencil or ppm_stencil.
• v_fs: A vector of functions to get the wind velocity at a given place and time. The function signature should be v_fs(i, j, k, t).
• Δ_fs: A vector of functions to get the grid spacing at a given place and time. The function signature should be Δ_fs(i, j, k, t).
• Δt: The time step size, which is assumed to be fixed.
• bc_type: The boundary condition type, e.g. ZeroGradBC().

Compute air mass for each layer from pressure edges and grid area.

Compute PBL mixing parameters based on domain levels.

Extract domain-specific pressure edges from the pre-computed 73-level GEOS-FP grid.

get_vf(domain, varname, data_f)

Return a function that gets the wind velocity at a given place and time for the given varname. data_f should be a function that takes a time and three spatial coordinates and returns the value of the wind speed in the direction indicated by varname.

get_Δ(domain, tff, pvaridx)

Return a function that gets the grid spacing at a given place and time for the given varname.

Apply full PBL mixing to tracer concentrations (GEOS-Chem TURBDAY algorithm).

pbl_obs_function(mtk_sys, coord_args, v, T)

Create an observed function for extracting data from the ModelingToolkit system.

resolve_species_bc(bc, x, species_vars)

Helper function to resolve species names to indices and create a SpeciesConstantBCArray. This is used by AdvectionOperator when species information is available.

Return the left and right stencil size of the first-order upwind stencil.

upwind1_stencil(ϕ, U, Δt, Δz; p)

First-order upwind advection in 1-D: https://en.wikipedia.org/wiki/Upwind_scheme.

• ϕ is the scalar field at the current time step, it should be a vector of length 3 (1 cell on the left, the central cell, and 1 cell on the right).
• U is the velocity at both edges of the central grid cell, it should be a vector of length 2.
• Δt is the length of the time step.
• Δz is the grid spacing.

The output will be time derivative of the central index (i.e. index 2) of the ϕ vector (i.e. dϕ/dt).

Δt and p are not used, but are function arguments for consistency with other operators.

Get a value from the x-direction velocity field.

Get a value from the y-direction velocity field.

Get a value from the z-direction velocity field.

function to get grid deltas.