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.

Zero gradient boundary conditions.

An array with zero gradient boundary conditions.

GaussianDispersion()

Return a ModelingToolkit.ODESystem 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).

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, DifferentialEquations

t0 = DateTime(2022, 5, 1)
t1 = DateTime(2022, 5, 2)

dom = DomainInfo(t0, t1; levrange=1:72,
                 lonrange = deg2rad(-125):deg2rad(5):deg2rad(-70),
                 latrange = deg2rad(25):deg2rad(4):deg2rad(54))

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

sys  = convert(ODESystem, mdl)

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

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

prob = ODEProblem(sys, u0, (datetime2unix(t0), datetime2unix(t1)), p)
sol = solve(prob, Tsit5();)

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

Puff(
    di::EarthSciMLBase.DomainInfo;
    buffer_cells,
    name
) -> ModelingToolkit.ODESystem

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)

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.

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.

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.