API Reference

surface_fluxes(
    param_set::APS,
    T_int,
    q_tot_int,
    q_liq_int,
    q_ice_int,
    ρ_int,
    T_sfc_guess,
    q_vap_sfc_guess,
    Φ_sfc,
    Δz,
    d,
    u_int = (0, 0),
    u_sfc = (0, 0),
    roughness_inputs = nothing,
    config = default_surface_flux_config(eltype(param_set)),
    scheme::SolverScheme = PointValueScheme(),
    solver_opts = nothing,
    flux_specs = nothing,
    update_T_sfc = nothing,
    update_q_vap_sfc = nothing,
)

Core entry point for calculating surface fluxes using Monin-Obukhov Similarity Theory (MOST).

Functionality

Calculates sensible heat flux, latent heat flux, momentum flux (stress), and friction velocity.

Can operate in four modes depending on inputs:

1. Prescribed Coefficients: If Cd and Ch are provided in flux_specs, fluxes are computed directly.
2. Fully Prescribed Fluxes: If shf, lhf, and ustar are provided, they are validated and the fluxes are returned.
3. Prescribed Heat and Drag: If shf, lhf, and Cd are provided, ustar is derived from Cd and wind speed.
4. Iterative Solver: Otherwise, iterates to find the Obukhov stability parameter ζ. Optional functions can be provided to calculate the skin temperature and skin humidity during the iteration.

Arguments

• param_set: SurfaceFluxes parameters (containing thermodynamics and universal function params).
• T_int: Interior (air) temperature [K] at height z.
• q_tot_int: Interior total specific humidity [kg/kg].
• q_liq_int, q_ice_int: Interior liquid/ice specific humidity [kg/kg].
• ρ_int: Interior air density [kg/m^3].
• T_sfc_guess: Initial guess for surface temperature [K], updated via callback if provided.
• q_vap_sfc_guess: Initial guess for surface vapor specific humidity [kg/kg], updated via callback if provided.
• Φ_sfc: Surface geopotential [m^2/s^2].
• Δz: Height of the reference (interior) level [m], measured from the surface under ReferenceAboveSurface or from the apparent sink d + z0m under ReferenceAboveApparentSink (set in config).
• d: Displacement height [m]. The Monin-Obukhov profiles span the effective height Δz - d above d, where the surface state applies (see surface_geopotential).
• u_int: Tuple of interior wind components (u, v) [m/s].
• u_sfc: Tuple of surface wind components (u, v) [m/s]. (Usually (0, 0)).
• roughness_inputs: Optional container of parameters (e.g., the plant area index PAI and canopy height h) that are passed directly to the specific roughness model (e.g., RaupachRoughnessParams).
• config: SurfaceFluxConfig struct containing:
  • roughness: Model for roughness lengths (e.g., ConstantRoughnessParams, COARE3RoughnessParams). Note: This package currently assumes the roughness length for heat (z0h) is equal to the roughness length for scalars (z0s).
  • gustiness: Model for gustiness (e.g., ConstantGustinessSpec).
  • moisture_model: DryModel or MoistModel.
  • rsl_model: Roughness sublayer model (e.g., NoRoughnessSubLayer, ExponentialRSL).
  • stability_cap: Cap on the stability parameter in stable conditions (e.g., NoStabilityCap, MaxHeatFluxStabilityCap).
  • reference_level: Convention for Δz (ReferenceAboveSurface or ReferenceAboveApparentSink).
• scheme: Discretization scheme (PointValueScheme or LayerAverageScheme).
• solver_opts: Options for the root solver (maxiter, tol, rtol, forced_fixed_iters).
• flux_specs: Optional FluxSpecs to prescribe specific constraints (e.g., ustar, shf, Cd).
• update_T_sfc: Optional callback update_T_sfc(ζ, param_set, thermo_params, inputs, scheme, u_star, z0m, z0h) that returns the surface temperature [K] during iteration.
• update_q_vap_sfc: Optional callback update_q_vap_sfc(ζ, param_set, thermo_params, inputs, scheme, T_sfc, u_star, z0m, z0h) that returns the surface vapor specific humidity [kg/kg] during iteration.

Returns

A SurfaceFluxConditions struct containing:

• shf: Sensible Heat Flux [W/m^2].
• lhf: Latent Heat Flux [W/m^2].
• evaporation: Evaporation rate [kg/m^2/s].
• ustar: Friction velocity [m/s].
• ρτxz, ρτyz: Momentum flux components (stress) [N/m^2].
• ζ: Stability parameter (z-d)/L [-].
• Cd: Drag coefficient [-].
• g_h: Heat conductance [m/s].
• T_sfc, q_vap_sfc: Final iterated surface temperature [K] and vapor specific humidity [kg/kg].
• L_MO: Monin-Obukhov length [m].
• L_eff, ζ_eff: Effective Obukhov length [m] and stability parameter [-] at which the exchange coefficients were evaluated (equal to L_MO and ζ unless a stability cap is active).
• converged: Convergence status; false, with all other fields NaN, when the reference level lies at or below a roughness length (see reference_height_valid and check_reference_height).

surface_fluxes(param_set, inputs, scheme=PointValueScheme(), solver_opts=nothing)

Dispatch to the appropriate solver mode based on the availability of inputs (coefficients, fluxes, or state).

SurfaceFluxConditions{FT}

Surface flux conditions returned by surface_fluxes.

All floating-point fields share the type FT, obtained by promoting the inputs. Momentum-flux components are the kinematic stress times density, i.e. ρτ = ρ u_* u_*, with units [kg/(m·s²)] = [N/m²].

Fields

• shf: Sensible heat flux [W/m²].
• lhf: Latent heat flux [W/m²].
• evaporation: Evaporation rate [kg/(m²·s)].
• ρτxz: Momentum flux, eastward component [kg/(m·s²)].
• ρτyz: Momentum flux, northward component [kg/(m·s²)].
• ustar: Friction velocity [m/s].
• ζ: Monin-Obukhov stability parameter (z - d)/L [-].
• Cd: Momentum exchange (drag) coefficient [-].
• g_h: Heat conductance Ch * U_eff [m/s].
• T_sfc: Surface temperature [K].
• q_vap_sfc: Surface air vapor specific humidity [kg/kg].
• L_MO: Monin-Obukhov length [m].
• L_eff: Effective Obukhov length for profile recovery, Δz_eff / min(ζ, ζ_cap) [m]. It equals L_MO unless a stability cap (see MaxHeatFluxStabilityCap) is active, in which case the exchange coefficients and similarity scales were evaluated at the capped stability parameter. Pass L_eff (not L_MO) to compute_profile_value to recover profiles consistent with the fluxes.
• ζ_eff: Stability parameter at which the exchange coefficients and similarity scales were evaluated, min(ζ, ζ_cap) = Δz_eff / L_eff [-]. It equals ζ unless a stability cap is active.
• converged: Solver convergence status. It is false when the reference level lies at or below a roughness length (see reference_height_valid), in which case all other fields are NaN.

The positional constructor accepts the fields in this order, with or without ζ_eff and L_eff (without ζ_eff, it is derived as ζ L_MO / L_eff; without L_eff, L_eff = L_MO).

SurfaceFluxConfig

Configuration for surface flux calculation components.

Fields

• roughness: The roughness length parameterization to use (e.g., ConstantRoughnessParams).
• gustiness: The gustiness parameterization to use (e.g., ConstantGustinessSpec).
• moisture_model: The moisture model (e.g., MoistModel or DryModel).
• rsl_model: Roughness sublayer correction model (e.g., ExponentialRSL). Defaults to NoRoughnessSubLayer (standard MOST, no RSL correction).
• stability_cap: Cap on the stability parameter in stable conditions (e.g., MaxHeatFluxStabilityCap). Defaults to NoStabilityCap (standard MOST).
• reference_level: Convention for the reference height Δz of the inputs, ReferenceAboveSurface (the default) or ReferenceAboveApparentSink.

FluxSpecs{FT}(; shf = nothing, lhf = nothing, ustar = nothing, Cd = nothing, Ch = nothing)
FluxSpecs(; shf = nothing, lhf = nothing, ustar = nothing, Cd = nothing, Ch = nothing)

Container for prescribed surface flux boundary conditions. Each field is nothing or a Real, including dual numbers for differentiation with respect to a prescribed value; surface_fluxes converts the values to the floating-point type of the state. The untyped constructor takes FT from the values (Float64 when none is given).

Fields

• shf: Sensible Heat Flux [W/m^2].
• lhf: Latent Heat Flux [W/m^2].
• ustar: Friction velocity [m/s].
• Cd: Momentum exchange coefficient.
• Ch: Heat exchange coefficient.

SolverOptions{FT}

Options for the Monin-Obukhov similarity theory solver.

Fields

• tol: Absolute tolerance on the stability parameter: the converged flag requires the final bracket width, or the step from the last iterate to the final regula falsi interpolant, to satisfy it, and in tolerance-checked mode it also bounds the step between iterates for the early exit.
• rtol: Relative tolerance on the stability parameter, used analogously to tol.
• maxiter: Number of bracket-refinement iterations. The ζ-solve performs 5 + maxiter residual evaluations in total (branch detection + bracketing probes + refinement); see the internal solve_stability_param.
• forced_fixed_iters: If true (default), disables the early tolerance exit and runs exactly maxiter refinement iterations (via RootSolvers.NoTolerance), so every point performs identical work (uniform control flow on GPUs). The converged flag is still evaluated from the final bracket (width and final step) and the tolerances.

SolverScheme

Abstract type for surface flux solver schemes.

PointValueScheme <: SolverScheme

Standard finite difference scheme using point values.

LayerAverageScheme <: SolverScheme

Finite volume approximation scheme following Nishizawa & Kitamura (2018).

Supported with BusingerParams and GryanikParams. Not supported with GrachevParams (use PointValueScheme instead).

compute_profile_value(param_set, L_MO, z0, Δz_eff, scale, val_sfc, transport, scheme, rsl_model)

Compute the (nondimensional) value of a variable (momentum or scalar) at effective aerodynamic height Δz_eff (height above surface minus displacement height).

Arguments

• param_set: Parameter set.
• L_MO: Monin-Obukhov length [m].
• z0: Roughness length [m].
• Δz_eff: Effective aerodynamic height z - d [m].
• scale: Similarity scale (ustar, thetastar, etc.).
• val_sfc: Surface value of the variable.
• transport: Transport type (MomentumTransport or HeatTransport).
• scheme: Discretization scheme (default: PointValueScheme()).
• rsl_model: Roughness sublayer model (default: NoRoughnessSubLayer). Use the same model as in the flux calculation for consistent profiles.

Formula:

X(Δz_eff) = (scale / κ) * F̂_z + val_sfc

where F̂_z = F_z + P is the dimensionless profile at height Δz_eff, including the roughness sublayer correction P (see rsl_corrected_profile).

screen_level_values(param_set, sc, inputs, z_screen, z_anemometer, scheme = PointValueScheme())

Return the NamedTuple (; T, q, u) of the air temperature [K] and vapor specific humidity [kg/kg] at the height z_screen [m] above the apparent sink for heat d + z0h, and the wind speed [m/s] at the height z_anemometer [m] above the apparent sink for momentum d + z0m, reconstructed from the Monin-Obukhov profiles of the solve that returned sc for inputs. The WMO screen and anemometer heights are 2 m and 10 m.

Between the surface and the reference level, a scalar that follows the heat profile takes the value $X_{sfc} + (X_{int} - X_{sfc}) r(z)$ with $r(z) = \widehat{F}_h(z) / \widehat{F}_h(Δz_{eff})$ (see dimensionless_profile_value), evaluated at the effective Obukhov length sc.L_eff, so that the profile reproduces the fluxes also under a stability cap. The temperature follows this relation in terms of the dry static energy, the variable the sensible heat flux is computed from, with the surface state at the displacement height (see surface_geopotential), and so includes the adiabatic change $g / c_{p,d}$ per meter between the screen and reference levels. The wind speed is $u_* \widehat{F}_m(z) / κ$, gustiness included, relative to the surface velocity u_sfc of the inputs.

The screen and anemometer values are point values of the profiles, while the reference profile $\widehat{F}_h(Δz_{eff})$ follows scheme. Under PointValueScheme, the profiles reach the interior state at the reference level, and the wind reaches the effective wind speed of the solve. Under LayerAverageScheme, the interior state is a layer average, which the point profile attains low in the layer, so point values in the upper part of the layer lie farther from the surface values than the interior state. Levels above the reference level take the values at the reference level. At or below the roughness length, the apparent sink at d + z0 above the surface, the dry static energy and humidity take the surface values and the wind vanishes, unless a roughness-sublayer model is configured: its correction keeps the profiles away from the surface values there.

Arguments

• param_set: Parameter set.
• sc: The SurfaceFluxConditions returned by surface_fluxes.
• inputs: The inputs container of that solve. See build_surface_flux_inputs.
• z_screen: Height of the screen level above the apparent sink for heat [m].
• z_anemometer: Height of the anemometer above the apparent sink for momentum [m].
• scheme: Discretization scheme of the solve (default: PointValueScheme).

Returns

A NamedTuple (; T, q, u):

• T: Air temperature at the screen height [K].
• q: Vapor specific humidity at the screen height [kg/kg].
• u: Wind speed at the anemometer height, relative to u_sfc [m/s].

dimensionless_profile_value(param_set, L_eff, z0, z, Δz_eff, transport, scheme, rsl_model)

Return the dimensionless Monin-Obukhov profile $\widehat{F}(z)$ of compute_profile_value at the height z above the displacement height, for the Obukhov length L_eff of a solve at the reference height Δz_eff. The height is clamped to the range from z0, where the profile is zero, to Δz_eff, so that the profile stays within the levels the solve connected.

Arguments

• param_set: Parameter set.
• L_eff: Effective Obukhov length of the solve, the field L_eff of SurfaceFluxConditions [m].
• z0: Roughness length of the transported quantity [m].
• z: Height above the displacement height [m].
• Δz_eff: Height of the reference level above the displacement height [m].
• transport: UF.MomentumTransport() or UF.HeatTransport().
• scheme: Discretization scheme (PointValueScheme or LayerAverageScheme).
• rsl_model: Roughness sublayer model of the solve (e.g., NoRoughnessSubLayer).

build_surface_flux_inputs(args...)

Centralized helper that normalizes user-facing specifications (winds, roughness, gustiness, flux constraints) into a NamedTuple containing inputs for surface flux calculations. Input types are preserved (no promotion), and users should ensure consistent input types for best performance.

Returns

A NamedTuple with the following fields:

Atmospheric and Surface State

• T_int: Interior air temperature [K]
• q_tot_int: Interior total specific humidity [kg/kg]
• q_liq_int: Interior liquid specific humidity [kg/kg]
• q_ice_int: Interior ice specific humidity [kg/kg]
• ρ_int: Interior air density [kg/m³]
• T_sfc_guess: Initial guess for surface temperature [K]. Can be nothing for default fallback.
• q_vap_sfc_guess: Initial guess for surface vapor specific humidity [kg/kg]. Can be nothing for default fallback.

Geometry

• Φ_sfc: Surface geopotential [m²/s²]
• Δz: Height of the reference level above the surface [m], under the convention reference_level
• d: Displacement height [m]
• reference_level: Convention for Δz, ReferenceAboveSurface or ReferenceAboveApparentSink; the solver converts the second to the first (see reference_above_surface)

Wind

• u_int: Horizontal wind components (u, v) at the interior level, as a tuple [m/s].
• u_sfc: Horizontal wind components (u, v) at the surface level, as a tuple [m/s].

Parameterizations

• roughness_model: Roughness parameterization, e.g. ConstantRoughnessParams.
• gustiness_model: Gustiness parameterization, e.g. ConstantGustinessSpec.
• moisture_model: Moisture model, MoistModel or DryModel.
• rsl_model: Roughness sublayer model, e.g. ExponentialRSL or NoRoughnessSubLayer.
• stability_cap: Stability cap specification, e.g. MaxHeatFluxStabilityCap or NoStabilityCap.
• ζ_cap: Numerical value of the stability cap, or nothing. It is nothing here and is set by the MOST solver from stability_cap (see with_stability_cap); functions that read the cap from the inputs compute it from stability_cap when it is nothing (see resolved_stability_cap).
• roughness_inputs: Optional inputs for roughness models.

Callbacks and Prescribed Values

• update_T_sfc: Optional callback to update surface temperature during iteration.
• update_q_vap_sfc: Optional callback to update surface vapor specific humidity during iteration.
• shf: Prescribed sensible heat flux from FluxSpecs; may be nothing [W/m²].
• lhf: Prescribed latent heat flux from FluxSpecs; may be nothing [W/m²].
• ustar: Prescribed friction velocity from FluxSpecs; may be nothing [m/s].
• Cd: Prescribed momentum exchange coefficient from FluxSpecs; may be nothing.
• Ch: Prescribed heat exchange coefficient from FluxSpecs; may be nothing.

shf = sensible_heat_flux(param_set, inputs, g_h, T_int, T_sfc, ρ_sfc, E)

Computes the sensible heat flux at the surface.

The sensible heat flux is given by

SHF = -ρ_sfc * g_h * ΔDSE + VSE_sfc * E

where ΔDSE = DSE_int - DSE_sfc is the difference in dry static energy between the interior and surface, g_h is the heat/moisture conductance, VSE_sfc is the dry static energy of water vapor at the surface temperature, and E is the evaporation rate. The second term, VSE_sfc * E, accounts for the vapor static energy VSE_sfc = cp_v * (T_sfc - T_0) + Φ_sfc (i.e., dry enthalpy cp_v * (T_sfc - T_0), or sensible heat, plus potential energy Φ_sfc) carried by evaporating water.

If inputs.shf is provided (not nothing), the function returns that value directly, allowing for prescribed sensible heat flux conditions. See the inputs container.

Arguments

• param_set: Parameter set.
• inputs: The inputs container. See build_surface_flux_inputs.
• g_h: Heat/moisture conductance [m/s].
• T_int: Interior temperature [K].
• T_sfc: Surface temperature [K].
• ρ_sfc: Surface air density [kg/m^3].
• E: Evaporation rate [kg/m^2/s]. Optional, default 0.

sensible_heat_flux(param_set, ζ, ustar, inputs, z0m, z0h, T_sfc, q_vap_sfc, ρ_sfc, scheme)

Computes the sensible heat flux given the Monin-Obukhov stability parameter ζ, friction velocity ustar, roughness lengths, and surface state. Useful for computing fluxes from variables available inside the solver loop.

Arguments

• param_set: Parameter set.
• ζ: Monin-Obukhov stability parameter.
• ustar: Friction velocity [m/s].
• inputs: The inputs container. See build_surface_flux_inputs.
• z0m: Momentum roughness length [m].
• z0h: Thermal roughness length [m].
• T_sfc: Surface temperature [K].
• q_vap_sfc: Surface vapor specific humidity [kg/kg].
• ρ_sfc: Surface air density [kg/m^3].
• scheme: Discretization scheme.

latent_heat_flux(param_set, inputs, E, model)

Computes the latent heat flux at the surface.

The latent heat flux is given by

LHF = LH_v0 * E

where LH_v0 is the latent heat of vaporization at the reference temperature and E is the evaporation rate.

If inputs.lhf is provided (not nothing), the function returns that value directly, allowing for prescribed latent heat flux conditions.

Arguments

• param_set: Parameter set.
• inputs: The inputs container. See build_surface_flux_inputs.
• E: Evaporation rate [kg/m^2/s].
• model: Moisture model (MoistModel or DryModel).

latent_heat_flux(param_set, ζ, ustar, inputs, z0m, z0h, q_vap_sfc, ρ_sfc, scheme)

Computes the latent heat flux given the Monin-Obukhov stability parameter ζ, friction velocity ustar, roughness lengths, and surface state. Calculates conductance and evaporation internally.

Arguments

• param_set: Parameter set.
• ζ: Monin-Obukhov stability parameter.
• ustar: Friction velocity [m/s].
• inputs: The inputs container. See build_surface_flux_inputs.
• z0m: Momentum roughness length [m].
• z0h: Thermal roughness length [m].
• q_vap_sfc: Surface vapor specific humidity [kg/kg].
• ρ_sfc: Surface air density [kg/m^3].
• scheme: Discretization scheme.

buoyancy_flux(param_set, shf, lhf, T_sfc, ρ_sfc, q_vap_sfc, q_liq_sfc, q_ice_sfc, model)

Computes the buoyancy flux at the surface, accounting for the presence of liquid and ice condensate.

The buoyancy flux B is defined as the vertical flux of virtual potential temperature θ_v. It is approximated by linearizing the density perturbations with respect to temperature and total water content:

B ≈ (g / ρ_sfc) * ( SHF / (cp_m * T_sfc) + (ε_vd - 1) * LHF / LH_v0 )

This form serves the prescribed-flux modes, where the fluxes are given and ζ is not. The MOST solve and its outputs use the exact form B = -u_*^3 ζ / (κ Δz_eff) of buoyancy_flux(param_set, ζ, ustar, inputs); the linearization differs from it by about a percent in strongly convective conditions.

Where:

• cp_m is the specific heat of moist air, calculated using q_tot_sfc, q_liq_sfc, and q_ice_sfc.
• ε_vd is the ratio of gas constants for water vapor and dry air.
• The term (ε_vd - 1) (approximately 0.61) represents the density effect of water vapor relative to dry air (the virtual temperature correction factor), ensuring the buoyancy flux accounts for the fact that moist air is lighter than dry air.

Arguments

• param_set: Parameter set.
• shf: Sensible heat flux [W/m²].
• lhf: Latent heat flux [W/m²].
• T_sfc: Surface temperature [K].
• ρ_sfc: Surface air density [kg/m³].
• q_vap_sfc: Surface water vapor specific humidity; default 0 [kg/kg].
• q_liq_sfc: Surface liquid water specific humidity; default 0 [kg/kg].
• q_ice_sfc: Surface ice specific humidity; default 0 [kg/kg].
• model: Moisture model (MoistModel or DryModel).

buoyancy_flux(param_set, ζ, ustar, inputs)

Computes the buoyancy flux given the Monin-Obukhov stability parameter ζ, friction velocity ustar, and geometric inputs via the inputs container.

The relationship is derived from the definition of the Obukhov length:

L = -u_*^3 / (κ * B)
ζ = Δz / L
=> B = -(u_*^3 * ζ) / (κ * Δz)

This is the buoyancy flux of the stability solve, consistent with its virtual potential temperatures; the prescribed-flux modes use the linearized form of buoyancy_flux(param_set, shf, lhf, T_sfc, ρ_sfc, ...) instead.

Arguments

• param_set: Parameter set.
• ζ: Monin-Obukhov stability parameter.
• ustar: Friction velocity [m/s].
• inputs: The inputs container. See build_surface_flux_inputs.

E = evaporation(param_set, inputs, g_h, q_vap_int, q_vap_sfc, ρ_sfc, model)

Computes the evaporation rate at the surface.

The evaporation rate is given by

E = -ρ_sfc * g_h * Δq_vap

where Δq_vap = q_vap_int - q_vap_sfc is the difference in vapor specific humidity between the interior and surface, g_h is the heat/moisture conductance (heat and moisture conductances must be equal for energetic consistency), and ρ_sfc is the surface air density. Here q_vap_int and q_vap_sfc are the vapor specific humidities (not total specific humidity) at the interior and surface, respectively.

If inputs.lhf is provided (not nothing), the function returns the evaporation rate computed from the prescribed latent heat flux: E = LHF / LH_v0, where LH_v0 is the latent heat of vaporization at the reference temperature.

Arguments

• param_set: Parameter set.
• inputs: The inputs container. See build_surface_flux_inputs.
• g_h: Heat conductance [m/s].
• q_vap_int: Interior vapor specific humidity [kg/kg].
• q_vap_sfc: Surface vapor specific humidity [kg/kg].
• ρ_sfc: Surface density [kg/m^3].
• model: Moisture model (MoistModel or DryModel).

evaporation(param_set, ζ, ustar, inputs, z0m, z0h, q_vap_sfc, ρ_sfc, scheme)

Computes the evaporation rate given the Monin-Obukhov stability parameter ζ, friction velocity ustar, roughness lengths, and surface state.

Arguments

• param_set: Parameter set.
• ζ: Monin-Obukhov stability parameter.
• ustar: Friction velocity [m/s].
• inputs: The inputs container. See build_surface_flux_inputs.
• z0m: Momentum roughness length [m].
• z0h: Thermal roughness length [m].
• q_vap_sfc: Surface vapor specific humidity [kg/kg].
• ρ_sfc: Surface air density [kg/m^3].
• scheme: Discretization scheme.

momentum_fluxes(Cd, inputs, ρ_sfc, gustiness)

Computes the momentum fluxes at the surface.

The momentum fluxes are calculated using the bulk aerodynamic formula:

ρτxz = -ρ_sfc * Cd * ΔU * Δu_x
ρτyz = -ρ_sfc * Cd * ΔU * Δu_y

where:

• Cd: Momentum exchange coefficient (drag coefficient)
• ΔU: Magnitude of the wind speed difference (including gustiness)
• Δu_x, Δu_y: Components of the wind speed difference
• ρ_sfc: Surface air density

Returns a tuple (ρτxz, ρτyz).

See the inputs container.

Arguments

• Cd: Drag coefficient.
• inputs: The inputs container. See build_surface_flux_inputs.
• ρ_sfc: Surface air density [kg/m^3].
• gustiness: Gustiness velocity scale [m/s].

state_bulk_richardson_number(param_set, inputs, T_sfc, ρ_sfc, ΔU, q_vap_sfc)

Computes the bulk Richardson number from the given state.

Arguments

• param_set: Parameter set.
• inputs: The inputs container. See build_surface_flux_inputs.
• T_sfc: Surface temperature [K].
• ρ_sfc: Surface air density [kg/m³].
• ΔU: Wind speed difference [m/s].
• q_vap_sfc: Surface vapor specific humidity [kg/kg]. Default: 0.

Returns the bulk Richardson number.

drag_coefficient(param_set, ζ, z0m, Δz_eff, scheme, rsl_model = NoRoughnessSubLayer())

Compute the drag coefficient Cd for momentum exchange.

Arguments

• param_set: Parameter set
• ζ: Stability parameter ζ = Δz_eff / L_MO
• z0m: Roughness length for momentum [m]
• Δz_eff: Effective aerodynamic height Δz - d [m]
• scheme: Surface flux solver scheme (default: PointValueScheme())
• rsl_model: Optional roughness sublayer model (default: NoRoughnessSubLayer).

Formula:

Cd = (κ / F̂_m)^2

where F̂_m = F_m + P_m is the RSL-corrected dimensionless velocity profile (see rsl_corrected_profile; F̂_m = F_m without RSL).

drag_coefficient(inputs, speed)

Compute the drag coefficient Cd from friction velocity (presumed to be in inputs.ustar) and effective wind speed (including any gustiness factors).

See the inputs container.

Arguments

• inputs: The inputs container. See build_surface_flux_inputs.
• speed: Effective wind speed [m/s].

heat_exchange_coefficient(param_set, ζ, z0m, z0h, Δz_eff, scheme, rsl_model = NoRoughnessSubLayer())

Compute the heat exchange coefficient Ch for scalar exchange.

Formula:

Ch = κ² / (F̂_m · F̂_h),

where F̂_m = F_m + P_m and F̂_h = F_h + P_h are the RSL-corrected dimensionless profiles for momentum and scalars respectively (see rsl_corrected_profile). For the finite-volume case, this corresponds to the formulation in Nishizawa & Kitamura (2018), Eqs. 21 & 22 (with Pr0 absorbed into Fh).

Arguments

• param_set: Parameter set
• ζ: Stability parameter ζ = Δz_eff / L_MO
• z0m: Roughness length for momentum [m]
• z0h: Roughness length for scalars (heat/moisture) [m]
• Δz_eff: Effective aerodynamic height Δz - d [m]
• scheme: Surface flux solver scheme (default: PointValueScheme())
• rsl_model: Optional roughness sublayer model (default: NoRoughnessSubLayer).

heat_conductance(param_set, ζ, ustar, inputs, z0m, z0h, scheme)

Compute the heat conductance g_h (speed * Ch), including any gustiness factor in the wind speed. Calculates windspeed and exchange coefficient internally from Monin-Obukhov variables; the gustiness uses the profile integrals of scheme (see gustiness_value). The exchange coefficient is evaluated at the stability parameter capped by any stability cap in inputs (see MaxHeatFluxStabilityCap and capped_stability).

Arguments

• param_set: Parameter set.
• ζ: Monin-Obukhov stability parameter.
• ustar: Friction velocity [m/s].
• inputs: The inputs container. See build_surface_flux_inputs.
• z0m: Momentum roughness length [m].
• z0h: Thermal roughness length [m].
• scheme: Discretization scheme.

compute_physical_scale_coeff(
    param_set::APS,
    Δz_eff,
    ζ,
    z0,
    transport,
    scheme::SolverScheme,
    rsl_model = NoRoughnessSubLayer(),
)

Compute the coefficient relating a bulk difference to its similarity scale.

Returns ϕ such that scale = Δvalue * ϕ; for example, u★ = ΔU * ϕ_m. It is given by

\[ϕ = \frac{κ}{\widehat{F}(Δz_{eff}, ζ, z_0)}\]

where κ is the von Kármán constant and F̂ = F + P is the RSL-corrected dimensionless profile (see rsl_corrected_profile; P = 0 when rsl_model is NoRoughnessSubLayer).

Arguments

• param_set: Parameter set.
• Δz_eff: Effective aerodynamic height Δz - d [m].
• ζ: Monin-Obukhov stability parameter [-].
• z0: Roughness length for the transported variable [m].
• transport: Transport type (MomentumTransport or HeatTransport).
• scheme: Discretization scheme (PointValueScheme or LayerAverageScheme).
• rsl_model: Optional roughness sublayer model (default: NoRoughnessSubLayer).

compute_ustar(param_set, ζ, z0, inputs, scheme, gustiness)

Return the friction velocity implied by the Monin-Obukhov solution.

If a friction velocity is prescribed via inputs.ustar (in the inputs container), it is returned directly; otherwise it is recomputed from the similarity coefficients, evaluated at the stability parameter capped by any stability cap in inputs (see MaxHeatFluxStabilityCap and capped_stability).

Arguments

• param_set: Parameter set.
• ζ: Monin-Obukhov stability parameter.
• z0: Momentum roughness length [m].
• inputs: The inputs container. See build_surface_flux_inputs.
• scheme: Discretization scheme.
• gustiness: Gustiness velocity scale [m/s].

compute_theta_star(param_set, ζ, z0h, inputs, scheme, T_sfc)

Return the potential temperature scale implied by the Monin-Obukhov solution, where z0h is the roughness length for heat.

Arguments

• param_set: Parameter set.
• ζ: Monin-Obukhov stability parameter.
• z0h: Thermal roughness length [m].
• inputs: The inputs container. See build_surface_flux_inputs.
• scheme: Discretization scheme.
• T_sfc: Surface temperature [K]. Optional; defaults to inputs.T_sfc_guess, falling back to the interior temperature inputs.T_int when the guess is nothing.

compute_q_star(param_set, ζ, z0h, inputs, scheme, q_vap_sfc)

Return the specific humidity scale implied by the current Monin-Obukhov solution, where z0h is the roughness length for scalars (assumed equal to heat).

Arguments

• param_set: Parameter set.
• ζ: Monin-Obukhov stability parameter.
• z0h: Thermal/scalar roughness length [m].
• inputs: The inputs container. See build_surface_flux_inputs.
• scheme: Discretization scheme.
• q_vap_sfc: Surface vapor specific humidity [kg/kg]. Optional; defaults to inputs.q_vap_sfc_guess, falling back to the interior total specific humidity inputs.q_tot_int when the guess is nothing.

surface_tke(param_set, Δz_eff, ustar, ζ)

Compute the surface-layer turbulent kinetic energy (TKE) (u_* ϕ)^2 following Tan et al. (2018).

Returns (u_* ϕ)^2 [m²/s²], where ϕ = sqrt(TKE)/u_* is the TKE-based velocity similarity function. In unstable conditions this is TKE = 3.75 u_*^2 + 0.2 w_*^2 + u_*^2 (-ζ)^{2/3}, and in stable/neutral conditions it reduces to 3.75 u_*^2. The convective (Deardorff) velocity scale w_* is computed from the mixed-layer height zi (a fixed parameter) and the buoyancy flux implied by ζ.

Arguments

• param_set: Parameter set.
• Δz_eff: Effective aerodynamic height Δz - d [m].
• ustar: Friction velocity [m/s].
• ζ: Monin-Obukhov stability parameter [-].

scalar_variance(param_set, scale, ζ)

Compute the scalar variance σ_s^2 = (scale * ϕ_σs)^2, using the temperature-variance similarity ϕ_σs = ϕ_σθ (Wyngaard et al. 1971; Tan et al. 2018).

Arguments

• param_set: Parameter set.
• scale: Similarity scale of the scalar (e.g., theta_star, q_star).
• ζ: Monin-Obukhov stability parameter.

theta_variance(param_set, inputs, shf, ustar, ζ, rho_sfc)

Computes potential temperature variance from sensible heat flux shf. Calculates θ_* = -shf / (ρ * c_p * u_*) and calls scalar_variance.

Arguments

• param_set: Parameter set.
• inputs: The inputs container. See build_surface_flux_inputs.
• shf: Sensible heat flux [W/m^2].
• ustar: Friction velocity [m/s].
• ζ: Monin-Obukhov stability parameter.
• rho_sfc: Surface density [kg/m^3].

obukhov_length(param_set, ustar, buoy_flux)

Computes the Monin-Obukhov length [m].

Returns zero if ustar is zero.

Arguments

• param_set: Parameter set.
• ustar: Friction velocity [m/s].
• buoy_flux: Surface buoyancy flux [m^2/s^3].

obukhov_stability_parameter(param_set, Δz_eff, ustar, buoy_flux)

Computes the Monin-Obukhov stability parameter ζ = Δz_eff / L_MO, where Δz_eff is the effective aerodynamic height ($z-d$).

Returns zero if ustar and hence L_MO are zero.

Arguments

• param_set: Parameter set.
• Δz_eff: Effective aerodynamic height [m].
• ustar: Friction velocity [m/s].
• buoy_flux: Surface buoyancy flux [m^2/s^3].

surface_density(param_set, T_int, ρ_int, T_sfc, Δz, q_tot_int=0, q_liq_int=0, q_ice_int=0, q_vap_sfc=nothing)
surface_density(param_set, inputs, T_sfc, q_vap_sfc)

Estimates the surface air density assuming hydrostatic balance between the interior and surface. It effectively extrapolates the interior pressure to the surface using the hydrostatic equation with an average virtual temperature, and then computes the surface density using the ideal gas law. The form with the inputs container extrapolates over the effective height Δz - d, from the reference level to the displacement height where the surface state applies (see surface_geopotential).

Arguments

• param_set: AbstractSurfaceFluxesParameters.
• T_int: Interior temperature [K].
• ρ_int: Interior density [kg/m^3].
• T_sfc: Surface temperature [K].
• Δz: Height difference [m].
• q_tot_int: Interior total specific humidity.
• q_liq_int: Interior liquid specific humidity.
• q_ice_int: Interior ice specific humidity.
• q_vap_sfc: Surface vapor specific humidity (optional, defaults to q_vap_int).

Returns ρ_sfc [kg/m^3].

surface_geopotential(param_set, inputs)

Compute the geopotential of the surface state, Φ_sfc + g * d [m²/s²]. The surface temperature and humidity apply at the apparent sink of the Monin-Obukhov profiles, which lies at the displacement height d above the surface (the roughness length z0h above it is neglected). Over a canopy, this is the level of the leaves that exchange heat with the air, so the dry static energy difference to the reference level, cp (T_int - T_sfc) + g (Δz - d), does not include the air column below the canopy.

Arguments

• param_set: Parameter set containing gravitational constant.
• inputs: The inputs container. See build_surface_flux_inputs.

surface_geopotential(inputs)

Return the geopotential of the ground, inputs.Φ_sfc [m²/s²]. Deprecated: the surface state applies at the displacement height, with the geopotential Φ_sfc + g d of surface_geopotential(param_set, inputs); a surface energy balance that uses this form with a displaced canopy is inconsistent with the sensible heat flux by g d / cp. Kept for callers of the one-argument form; to be removed in the next breaking release.

interior_geopotential(param_set, inputs)

Compute the geopotential at the interior (atmospheric) reference level.

Arguments

• param_set: Parameter set containing gravitational constant.
• inputs: The inputs container. See build_surface_flux_inputs.

Returns Φ_sfc + g * Δz [m²/s²], with Δz measured from the surface (see reference_above_surface).

effective_height(inputs)
effective_height(param_set, inputs)

Compute the effective aerodynamic height z_eff = Δz - d, the height of the reference level above the displacement height, which the Monin-Obukhov profiles span and over which the surface state at d (see surface_geopotential) is connected to the interior state. The one-argument form expects inputs under ReferenceAboveSurface; the two-argument form converts inputs under ReferenceAboveApparentSink first with reference_above_surface.

Arguments

• param_set: Parameter set (required when inputs may use ReferenceAboveApparentSink).
• inputs: The inputs container. See build_surface_flux_inputs.

Returns Δz - d [m].

ReferenceAboveSurface

The reference height Δz of the inputs is measured from the surface, and the Monin-Obukhov profiles span the effective height Δz - d above the displacement height d. This is the default convention.

ReferenceAboveApparentSink

The reference height Δz of the inputs is measured from the apparent sink for momentum, d + z0m above the surface. The solve converts it to the height Δz + d + z0m above the surface, so the reference level lies above the sink for any displacement height, and the profiles span Δz + z0m. Land models forced by reanalysis or by an atmosphere model that does not resolve the canopy use this convention: their forcing heights are defined relative to the surface the atmosphere feels, not to the ground below a canopy.

The conversion uses the roughness length before the solve, so the roughness model must be independent of the friction velocity (see depends_on_ustar and reference_above_surface).

reference_above_surface(param_set, inputs)

Return the inputs with the reference height Δz measured from the surface. Under ReferenceAboveSurface, the inputs are returned as they are. Under ReferenceAboveApparentSink, Δz is measured from the apparent sink for momentum and becomes Δz + d + z0m, with the roughness length z0m of a roughness model that is independent of the friction velocity.

Arguments

• param_set: Parameter set.
• inputs: The inputs container. See build_surface_flux_inputs.

interior_vapor_specific_humidity(inputs)

Return the vapor specific humidity of the interior air [kg/kg], the total specific humidity q_tot_int less the condensate q_liq_int + q_ice_int.

Arguments

• inputs: The inputs container. See build_surface_flux_inputs.

reference_height_valid(inputs, z0m, z0h = z0m)

Whether the reference level lies above both roughness lengths, Δz - d > max(z0m, z0h), so that the Monin-Obukhov profiles of momentum and of scalars between the surface and the reference level are defined. The scalar roughness length matters when it exceeds the momentum one, as the COARE 3.0 model gives at low friction velocities. Δz is the height above the surface: inputs under ReferenceAboveApparentSink are converted first with reference_above_surface, as surface_fluxes does. surface_fluxes returns NaN fluxes with converged = false for inputs that fail this test, since the solve cannot throw inside a GPU kernel; check_reference_height raises the corresponding error on the host.

Arguments

• inputs: The inputs container. See build_surface_flux_inputs.
• z0m: Momentum roughness length [m].
• z0h: Scalar roughness length [m]; z0m by default.

check_reference_height(Δz, d, z0m, z0h = z0m)

Throw an ArgumentError unless the reference level lies above both roughness lengths, Δz - d > max(z0m, z0h) (see reference_height_valid). For a host-side check of a model's configuration before fluxes are computed in kernels.

Arguments

• Δz: Height of the reference level above the surface [m].
• d: Displacement height [m].
• z0m: Momentum roughness length [m].
• z0h: Scalar roughness length [m]; z0m by default.

invalidate_unless(sc::SurfaceFluxConditions, valid::Bool)

Return sc when valid, and otherwise a copy with every floating-point field set to NaN and converged = false. Written with ifelse so that it runs in GPU kernels, where the solve cannot throw; see reference_height_valid.

ConstantRoughnessParams{FT} <: AbstractRoughnessParams

Roughness lengths fixed to constant values.

Fields

• z0m: Momentum roughness length [m].
• z0s: Scalar roughness length [m]. Used for both heat (z0h) and humidity (z0q).

The keyword defaults shown here (z0m = 2e-4 m, z0s = 2e-5 m) are also the roughness lengths used by surface_fluxes when config is omitted, via default_surface_flux_config. When loading via ClimaParams, the values are read from the TOML file. Most applications pass z0m/z0s explicitly or load them from ClimaParams.

COARE3RoughnessParams{FT} <: AbstractRoughnessParams

COARE 3.0 roughness parameterization (Fairall et al. 2003).

References

• Fairall, C. W., Bradley, E. F., Hare, J. E., Grachev, A. A., & Edson, J. B. (2003). Bulk parameterization of air–sea fluxes: Updates and verification for the COARE algorithm. Journal of Climate, 16, 571–591. DOI: 10.1175/1520-0442(2003)016<0571:BPOASF>2.0.CO;2

The default values specified here are used when constructing the struct manually. When loading via ClimaParams, these values are overwritten by the parameters in the TOML file.

RaupachRoughnessParams <: AbstractRoughnessParams

Raupach (1994) canopy roughness model: the momentum roughness length z0m and the zero-plane displacement height d as functions of the canopy height h and the area index of the canopy (see momentum_roughness and displacement_height). The roughness inputs are the canopy height roughness_inputs.h and the plant area index Λ = roughness_inputs.PAI, the single-sided area of all canopy elements (leaves, living or dead, stems, and branches) per unit ground area, which Raupach (1994) calls the canopy area index: the sum of the leaf and stem area indices, LAI + SAI. The field name LAI is accepted as a deprecated alias for PAI.

Raupach (1994) writes the drag partition in terms of the frontal area index λ, the frontal area of the canopy elements facing the mean wind per unit ground area, with Λ = 2λ for isotropically oriented elements (see frontal_area_index and canopy_area_index).

Fields

• C_R: Drag coefficient of an isolated roughness element (0.3).
• C_S: Drag coefficient of the substrate at height h (0.003).
• c_d1: Constant of the displacement height expression (7.5).
• stanton_number: Ratio z0s / z0m of the scalar to the momentum roughness length; exp(-kB⁻¹) for an excess resistance kB⁻¹ = ln(z0m / z0s) (0.1, so that kB⁻¹ ≈ 2.3).
• frontal_area_ratio: Frontal area index per unit plant area index, λ = frontal_area_ratio * Λ; 0.5 for isotropically oriented elements (Raupach 1994). It enters the drag partition (Eq. 7) only; the displacement height (Eq. 8) is an empirical fit in Λ.
• λ_min: Floor on the frontal area index (0, so that PAI is used as given, following Raupach 1994). A positive floor stands in for stems and branches when the input counts leaves only, or vanishes for a canopy of nonzero height, so that such a canopy stays aerodynamically rough; with a plant area index that includes the stem area, it is unnecessary.
• ustar_Uh_max: Sheltering limit of u★ / U(h) (0.3, Raupach 1994, Eq. 7).
• c_w: Ratio (z_w - d) / (h - d) of the heights of the roughness-sublayer top z_w and the canopy top above the displacement height (2, Raupach 1994). It sets the roughness-sublayer influence function at the canopy top, Ψ_h = ln c_w - 1 + 1 / c_w (Eq. 5), which is 0.193 for c_w = 2.

References

• Raupach, M. R. (1994). Simplified expressions for vegetation roughness length and zero-plane displacement as functions of canopy height and area index. Boundary-Layer Meteorology, 71, 211–216. DOI: 10.1007/BF00709229

The default values specified here are used when constructing the struct manually. When loading via ClimaParams, all fields are read from the TOML file: stanton_number and the raupach_* parameters.

momentum_roughness(spec::COARE3RoughnessParams, u★, sfc_param_set, roughness_inputs)

Calculate momentum roughness length using the COARE 3.0 algorithm (Fairall et al. 2003).

Formulation

The momentum roughness length z0m is parameterized as the sum of a smooth flow limit (Smith 1988) and a rough flow limit (Charnock 1955):

\[z_{0m} = z_{0m,smooth} + z_{0m,rough}\]

• Smooth flow: Dominated by viscous limit, proportional to ν / u★.
• Rough flow: Dominated by wind stress, proportional to α * u★^2 / g.

The Charnock parameter α varies with the 10-m wind speed and is interpolated linearly between lower and upper bounds defined in spec.

Dependencies

• u★: Friction velocity [m/s]
• kinematic_visc: Kinematic viscosity of air [m^2/s], from spec
• grav: Gravitational acceleration [m/s^2], from sfc_param_set
• mag_u_10: 10m wind speed [m/s]. (internally recovered from u★ assuming neutral log profile)

References

• Fairall, C. W., Bradley, E. F., Hare, J. E., Grachev, A. A., & Edson, J. B. (2003). Bulk parameterization of air–sea fluxes: Updates and verification for the COARE algorithm. Journal of Climate, 16, 571–591. DOI: 10.1175/1520-0442(2003)016<0571:BPOASF>2.0.CO;2
• Smith, S. D. (1988). Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature. Journal of Geophysical Research: Oceans, 93, 15467–15472. DOI: 10.1029/JC093iC12p15467
• Charnock, H. (1955). Wind stress on a water surface. Quarterly Journal of the Royal Meteorological Society, 81, 639–640. DOI: 10.1002/qj.49708135027

momentum_roughness(spec::RaupachRoughnessParams, u★, sfc_param_set, roughness_inputs)

Momentum roughness length [m] of the Raupach (1994) canopy roughness model, h times raupach_roughness_fraction at the plant area index roughness_inputs.PAI, and at least the fixed roughness length z0m_fixed of the parameter set (which also covers a vanishing canopy height).

Formulation

The model partitions the surface drag between the substrate (soil) and the roughness elements (plants), which gives the friction velocity ratio at the canopy top u★ / U(h) = min((C_S + C_R λ)^(1/2), (u★ / U(h))_max) for the frontal area index λ (Eq. 7); the cap is the sheltering limit, beyond which z0m / h decreases with λ (for λ > 0.29, i.e., PAI > 0.58 with frontal_area_ratio = 0.5). The roughness length follows from the wind profile at the canopy top with the roughness-sublayer influence function Ψ_h (Eq. 4), set by the roughness-sublayer depth ratio c_w (Eq. 5).

Dependencies

• roughness_inputs.PAI: Plant area index Λ = LAI + SAI (leaves plus stems) [m^2/m^2], from which the frontal area index is λ = max(frontal_area_ratio * Λ, λ_min). Λ enters the displacement height (Eq. 8) directly and λ the drag partition (Eq. 7).
• roughness_inputs.h: Canopy height [m].
• spec: Coefficients, see RaupachRoughnessParams.
• z0m_fixed and von_karman_const from sfc_param_set.

References

• Raupach, M. R. (1994). Simplified expressions for vegetation roughness length and zero-plane displacement as functions of canopy height and area index. Boundary-Layer Meteorology, 71, 211–216. DOI: 10.1007/BF00709229

scalar_roughness(spec::COARE3RoughnessParams, u★, sfc_param_set, roughness_inputs)

Calculate scalar roughness length using the COARE 3.0 algorithm (Fairall et al. 2003).

Formulation

The scalar roughness length z0s is parameterized as an empirical fit to COARE and HEXOS data. It limits z0s to a smooth flow limit (1.1e-4 m) and decreases for rough flow following a power law of the roughness Reynolds number (Re_star):

\[z_{0s} = \min(1.1 \times 10^{-4}, 5.5 \times 10^{-5} \cdot R_{e*}^{-0.6})\]

where $R_{e*} = z_{0m} u_* / \nu$.

Dependencies

• u★: Friction velocity [m/s]
• kinematic_visc: Kinematic viscosity of air [m^2/s]
• z0m: Momentum roughness length (calculated internally)

References

• Fairall, C. W., Bradley, E. F., Hare, J. E., Grachev, A. A., & Edson, J. B. (2003). Bulk parameterization of air–sea fluxes: Updates and verification for the COARE algorithm. Journal of Climate, 16, 571–591. DOI: 10.1175/1520-0442(2003)016<0571:BPOASF>2.0.CO;2

scalar_roughness(spec::RaupachRoughnessParams, u★, sfc_param_set, roughness_inputs)

Calculate scalar roughness length from the momentum roughness length using a fixed Stanton number.

Formulation

Scaled from the momentum roughness length using a fixed Stanton number:

\[z_{0s} = z_{0m} \cdot St\]

where $St$ is spec.stanton_number.

Dependencies

• spec.stanton_number: Stanton number specific to the canopy/surface type.

displacement_height(spec::RaupachRoughnessParams, roughness_inputs)

Zero-plane displacement height [m] of the canopy, h times raupach_displacement_fraction at the plant area index roughness_inputs.PAI. Inside surface_fluxes, the displacement height is the input d, which sets the effective reference height Δz - d; callers compute it here from the same canopy inputs as the roughness length.

frontal_area_index(spec::RaupachRoughnessParams, plant_area_index)

Frontal area index λ = max(frontal_area_ratio * plant_area_index, λ_min) [m^2/m^2] of the canopy, which sets the drag partition of the Raupach (1994) roughness model (Eq. 7).

canopy_area_index(spec::RaupachRoughnessParams, plant_area_index)

Canopy area index Λ [m^2/m^2] of the Raupach (1994) displacement height (Eq. 8): the plant area index, raised to λ_min / frontal_area_ratio where the floor on the frontal_area_index is active, so that both indices describe the same canopy.

raupach_displacement_fraction(spec::RaupachRoughnessParams, plant_area_index)

Zero-plane displacement height as a fraction of the canopy height for a plant area index (Raupach 1994, Eq. 8), with the canopy_area_index Λ:

\[d / h = 1 - \frac{1 - \exp(-\sqrt{c_{d1} Λ})}{\sqrt{c_{d1} Λ}}\]

The fraction tends to zero as Λ → 0.

raupach_roughness_fraction(spec::RaupachRoughnessParams, κ, plant_area_index)

Momentum roughness length as a fraction of the canopy height for a plant area index and von Kármán constant κ (Raupach 1994, Eqs. 4, 5, 7, and 8):

\[z_{0m} / h = (1 - d / h) \exp(-κ U_h / u_* - Ψ_h), \quad u_* / U_h = \min(\sqrt{C_S + C_R λ}, (u_* / U_h)_{max}),\]

with the frontal_area_index λ, d / h from raupach_displacement_fraction, and the roughness-sublayer influence function Ψ_h = ln c_w - 1 + 1 / c_w, the departure of the wind profile immediately above the canopy from the logarithmic law. The cap on u_* / U_h is the sheltering limit, beyond which z0m / h decreases with λ (for λ > 0.29 with the default coefficients).

Arguments

• spec: Coefficients, see RaupachRoughnessParams.
• κ: Von Kármán constant [-].
• plant_area_index: Plant area index, the sum of the leaf and stem area indices [m^2/m^2].

ConstantGustinessSpec{TG <: Real}

A gustiness model where the gustiness velocity is a constant value.

Fields

• value: The constant gustiness velocity [m/s].

DeardorffGustinessSpec

A gustiness model based on Deardorff (1970) scaling with convective velocity scale $w_*$. The gustiness velocity is computed as: $U_{gust} = \beta w_*$ where $w_* = (B z_i)^{1/3}$.

FlooredDeardorffGustinessSpec{FT <: Real}

The larger of a constant minimum wind speed u_min [m/s] and the convective gustiness $β w_*$ of DeardorffGustinessSpec, with the convective velocity scale $w_* = (B z_i)^{1/3}$ of Deardorff (1970) for a positive surface buoyancy flux $B$, the boundary layer depth $z_i$ (gustiness_zi) and coefficient $β$ (gustiness_coeff). The convective part vanishes in stable conditions, where the floor applies; FlooredDeardorffGustinessSpec(zero(FT)) is the pure convective gustiness.

Within the stability solve, the convective part is evaluated in closed form from the surface and atmospheric state (see free_convection_wind_speed): at a stability parameter $ζ$, the bulk relations $u_* = κ U / F_m$ and $θ_{v*} = κ Δθ_v / F_h$ make the buoyancy flux $B = (g/θ_v) u_* θ_{v*}$ linear in the effective wind speed $U$, so that $U = β w_*$ has the solution

\[U^2 = β^3 κ^2 \frac{g}{θ_v} z_i \frac{Δθ_v}{F_m(ζ) F_h(ζ)}.\]

The effective wind speed is the fixed point of $U = \max(|Δu|, u_{min}, β w_*(U))$, and the gustiness is independent of $u_*$ (see depends_on_ustar), so the friction velocity follows from $ζ$ in closed form and the free-convection limit is well posed at every $ζ$. DeardorffGustinessSpec evaluates the gustiness from the buoyancy flux implied by $ζ$ and the current $u_*$; at fixed $ζ$ that gustiness is proportional to $u_*$, and a consistent $u_*$ exists only up to the free-convection limit.

Fields

• u_min: Minimum wind speed [m/s].

Examples

gustiness = FlooredDeardorffGustinessSpec(0.5)  # floor of 0.5 m/s
config = SurfaceFluxConfig(ConstantRoughnessParams(0.01, 0.001), gustiness)

References

• Deardorff, J. W. (1970). Convective velocity and temperature scales for the unstable planetary boundary layer and for Rayleigh convection. J. Atmos. Sci., 27, 1211-1213.
• Beljaars, A. C. M. (1995). The parametrization of surface fluxes in large-scale models under free convection. Q. J. R. Meteorol. Soc., 121, 255-270.

MoistModel

Indicates that moisture effects (latent heat, virtual temperature) should be included in the flux calculations.

DryModel

Indicates that moisture effects should be ignored (sensible heat and momentum only).

gustiness_value(spec, param_set, buoyancy_flux)

Returns the gustiness velocity scale [m/s] based on the specification.

Arguments

• spec: The gustiness specification (e.g., ConstantGustinessSpec or DeardorffGustinessSpec).
• param_set: Parameter set containing constants and coefficients.
• buoyancy_flux: Surface buoyancy flux [m^2/s^3], required for Deardorff gustiness.

gustiness_value(::DeardorffGustinessSpec, param_set, buoyancy_flux)

Calculates the gustiness based on the Deardorff convective velocity scale.

Formulation

The gustiness $U_{gust}$ is parameterized as proportional to the Deardorff velocity $w_*$:

\[U_{gust} = C_{gust} \cdot w_*\]

where $w_* = (B \cdot z_i)^{1/3}$.

• $B$ is the surface buoyancy flux (buoyancy_flux).
• $z_i$ is the boundary layer height (assumed fixed togustiness_zi from parameters).
• $C_{gust}$ is a scaling coefficient (gustiness_coeff from parameters).

This formulation parametrizes the enhancement of surface fluxes due to boundary layer scale eddies in unstable conditions, particularly important in low-wind regimes (free convection limit).

References

• Deardorff, J. W. (1970). Convective velocity and temperature scales for the unstable planetary boundary layer and for Rayleigh convection. Journal of the Atmospheric Sciences, 27, 1211-1213. DOI: 10.1175/1520-0469(1970)027<1211:CVATSF>2.0.CO;2
• Beljaars, A. C. M. (1995). The parametrization of surface fluxes in large-scale models under free convection Quarterly Journal of the Royal Meteorological Society, 121, 255-270. DOI: 10.1002/qj.49712152203

gustiness_value(spec, param_set, ζ, ustar, inputs, scheme = PointValueScheme())

Return the gustiness velocity scale [m/s] from the solver variables ζ and ustar. ConstantGustinessSpec returns its value; FlooredDeardorffGustinessSpec evaluates its convective part in closed form with the profile integrals of scheme; the other models evaluate the buoyancy flux first (see depends_on_ustar).

gustiness_value(spec::FlooredDeardorffGustinessSpec, param_set, buoyancy_flux)

Return the larger of the floor spec.u_min and the convective gustiness $β w_*$ of DeardorffGustinessSpec for the surface buoyancy flux buoyancy_flux [m²/s³], zero for a non-positive buoyancy flux. The post-solve fluxes use this form, with the buoyancy flux implied by the converged ζ and friction velocity.

gustiness_value(
    spec::FlooredDeardorffGustinessSpec, param_set, ζ, ustar, inputs,
    scheme = PointValueScheme(),
)

Return the larger of the floor spec.u_min and the free-convection wind speed free_convection_wind_speed at the stability parameter ζ, with the profile integrals of scheme. The friction velocity ustar enters only through the roughness lengths.

minimum_wind_speed(spec::AbstractGustinessSpec, param_set)

Return the minimum effective wind speed [m/s] that a gustiness model imposes in all conditions, in the floating-point type of param_set: the value of a ConstantGustinessSpec, the floor u_min of a FlooredDeardorffGustinessSpec, and zero for models whose gustiness vanishes in stable conditions, such as DeardorffGustinessSpec. A model that folds this floor into the wind speed it passes to the solve (for example, a canopy model that attenuates the wind above the canopy to the ground below it) pairs the attenuated wind with without_floor.

without_floor(spec::AbstractGustinessSpec)

Return the gustiness model spec with its minimum wind speed set to zero: a ConstantGustinessSpec becomes a zero gustiness, a FlooredDeardorffGustinessSpec keeps its convective part only, and models with no floor are returned as they are. See minimum_wind_speed.

free_convection_wind_speed(param_set, ζ, ustar, inputs, scheme = PointValueScheme())

Return the effective wind speed $U$ [m/s] at which the convective gustiness $β w_*$ equals $U$ itself, for the Monin-Obukhov stability parameter ζ and the surface and atmospheric state in inputs:

\[U^2 = β^3 κ^2 \frac{g}{θ_v} z_i \frac{Δθ_v}{F_m(ζ) F_h(ζ)},\]

where $Δθ_v$ is the virtual potential temperature excess of the surface over the air, $θ_v$ the air value, $F_m$ and $F_h$ the dimensionless profile integrals of momentum and heat (so that $u_* = κ U / F_m$ and $θ_{v*} = κ Δθ_v / F_h$), and $β$, $z_i$, $κ$, $g$ the parameters gustiness_coeff, gustiness_zi, von_karman_const, and grav. The result follows from $w_*^3 = B z_i$ with $B = (g / θ_v) u_* θ_{v*}$ and $U = β w_*$. It is zero when the surface is not warmer than the air ($Δθ_v ≤ 0$). The profile integrals are evaluated with the discretization scheme of the solve (PointValueScheme or LayerAverageScheme), the roughness sublayer model, and the stability cap of inputs, with the roughness lengths of the roughness model at ustar. The surface temperature and humidity are the guesses T_sfc_guess and q_vap_sfc_guess of inputs. Without surface callbacks these are the values the stability solve uses in its bulk Richardson number, so the closed form is the exact fixed point of the solver's bulk relations at ζ. With callbacks, the solve advances the guesses between residual evaluations, so the gustiness lags the surface state of the bulk Richardson number by one evaluation, as the friction velocity does; the two agree once the surface state has converged.

virtual_pottemps(param_set, inputs, T_sfc, ρ_sfc, q_vap_sfc = 0)

Return the virtual potential temperatures (θ_v_sfc, θ_v_int) [K] of the surface and of the interior air. Called from state_bulk_richardson_number and free_convection_wind_speed, so both evaluate the same $Δθ_v$. The condensate concentration is taken to be the same at the surface and in the interior.

Arguments

• param_set: Parameter set.
• inputs: The inputs container. See build_surface_flux_inputs.
• T_sfc: Surface temperature [K].
• ρ_sfc: Surface air density [kg/m³].
• q_vap_sfc: Surface vapor specific humidity [kg/kg]. Default: 0.

depends_on_ustar(model)

Whether the roughness lengths of a roughness model, or the gustiness of a gustiness model, depend on the friction velocity u★. compute_ustar_and_roughness finds u★ by root-finding when either model does; otherwise, the roughness lengths, the gustiness, and u★ follow directly from the stability parameter. The models independent of u★ are ConstantRoughnessParams, RaupachRoughnessParams, ConstantGustinessSpec, and FlooredDeardorffGustinessSpec.

compute_ustar_and_roughness(param_set, ζ, inputs, scheme)

Computes friction velocity ustar and roughness lengths z0m, z0h for a given stability ζ.

• If inputs.ustar is prescribed, it is returned directly.
• If the roughness and gustiness models are independent of ustar (see depends_on_ustar), z0m, z0h, and ustar follow directly from ζ.
• Otherwise, three iterations of Brent's method on ustar ∈ [1e-4, 4] m/s find the ustar consistent with the roughness and gustiness models; the result lies in this bracket. If no consistent ustar lies in the bracket, the endpoint on the side of the root is returned: 4 m/s when the friction velocity implied by ζ and the gustiness it generates exceeds the bracket for every ustar, and 1e-4 m/s in calm conditions. With DeardorffGustinessSpec, the gustiness at fixed ζ is proportional to ustar, and no consistent ustar exists for ζ more unstable than the free-convection limit; the solve for ζ then settles where a consistent ustar exists.

NoRoughnessSubLayer <: AbstractRoughnessSubLayerModel

No roughness sublayer correction. Standard Monin-Obukhov similarity theory (MOST) is applied without modification. This is the default when no RSL model is specified.

LinearRSL{FT} <: AbstractRoughnessSubLayerModel
LinearRSL(; c_m = 0.4, c_h = 0.4, z_RSL = 10.0)
LinearRSL(FT; kwargs...)

Roughness sublayer model with a linear RSL factor,

μ(z) = 1 - c (1 - z / z_RSL)   for z < z_RSL,     μ(z) = 1   for z ≥ z_RSL,

where z is the height above the displacement height d, so that the dimensionless gradients are φ̂(z) = φ(z / L) μ(z). The factor increases linearly from 1 - c at the displacement height to 1 at the top of the RSL. It is the first-order (small-c) approximation of ExponentialRSL.

The corrected profile coincides with the MOST profile above the RSL, so the roughness length z0 and displacement height d are the apparent values, as obtained from standard canopy relations (e.g., z0 ≈ 0.1 h, d ≈ 0.67 h, or RaupachRoughnessParams) or by fitting MOST profiles above the RSL. Within the RSL, the wind speed and scalar differences are larger, and the exchange coefficients smaller, than those of the MOST profile extrapolated downward with the apparent z0 and d.

Fields

• c_m: RSL strength for momentum, μ(0) = 1 - c_m, with 0 ≤ c_m < 1 [-].
• c_h: RSL strength for scalars (heat, moisture), with 0 ≤ c_h < 1 [-].
• z_RSL: RSL depth above the displacement height, ≥ 0 [m].

The RSL top is typically at 2–3 canopy heights h above the ground (Garratt 1980; Raupach et al. 1991), i.e., z_RSL ≈ 1.3–2.3 h for d ≈ 0.67 h; Physick & Garratt (1995) use a depth of 50 z0. z_RSL is a fixed depth, to be scaled with the canopy height when the model is constructed.

FT is the floating-point type of the parameters (default Float64). In the flux computation, the parameters are converted to the floating-point type of the inputs, so Float32 computations remain in Float32 with the default constructors; LinearRSL(Float32; ...) constructs a model with Float32 parameters.

Examples

h = 30.0  # canopy height [m]
rsl = LinearRSL(c_m = 0.4, c_h = 0.4, z_RSL = 2h - 0.67h)
rsl32 = LinearRSL(Float32; z_RSL = 20.0)

References

• Garratt, J. R. (1980). Surface influence upon vertical profiles in the atmospheric near-surface layer. Quarterly Journal of the Royal Meteorological Society, 106, 803–819.
• Physick, W. L., & Garratt, J. R. (1995). Incorporation of a high-roughness lower boundary into a mesoscale model for studies of dry deposition over complex terrain. Boundary-Layer Meteorology, 74, 55–71.
• Raupach, M. R., Antonia, R. A., & Rajagopalan, S. (1991). Rough-wall turbulent boundary layers. Applied Mechanics Reviews, 44, 1–25.

ExponentialRSL{FT} <: AbstractRoughnessSubLayerModel
ExponentialRSL(; c_m = 0.7, c_h = 0.7, z_RSL = 10.0)
ExponentialRSL(FT; kwargs...)

Roughness sublayer model with an exponential RSL factor,

μ(z) = exp(-c (1 - z / z_RSL))   for z < z_RSL,     μ(z) = 1   for z ≥ z_RSL,

where z is the height above the displacement height d, so that the dimensionless gradients are φ̂(z) = φ(z / L) μ(z). The factor increases from exp(-c) at the displacement height to 1 at the top of the RSL. This is the form of Garratt (1980, 1983), used by Physick & Garratt (1995) with c = 0.7 (their 0.5 exp(0.7 z / z_RSL), as ln 2 ≈ 0.7) and shown in Bonan (2019, Fig. 6.8). Physick & Garratt (1995) also anchor the corrected profiles at the RSL top (their Eqs. 7 and 9), as done here.

The RSL factor μ is a prescribed function of height; the stability dependence of the corrected profiles enters through φ(z / L). As for LinearRSL, the corrected profile coincides with MOST above the RSL, so z0 and d are the apparent roughness length and displacement height.

Fields

• c_m: RSL exponent for momentum, μ(0) = exp(-c_m), with c_m ≥ 0 [-].
• c_h: RSL exponent for scalars (heat, moisture), with c_h ≥ 0 [-].
• z_RSL: RSL depth above the displacement height, ≥ 0 [m].

See LinearRSL for typical RSL depths and the floating-point type.

Examples

h = 30.0  # canopy height [m]
rsl = ExponentialRSL(c_m = 0.7, c_h = 0.7, z_RSL = 2h - 0.67h)

References

• Garratt, J. R. (1980). Surface influence upon vertical profiles in the atmospheric near-surface layer. Quarterly Journal of the Royal Meteorological Society, 106, 803–819.
• Physick, W. L., & Garratt, J. R. (1995). Incorporation of a high-roughness lower boundary into a mesoscale model for studies of dry deposition over complex terrain. Boundary-Layer Meteorology, 74, 55–71.
• Harman, I. N., & Finnigan, J. J. (2007). A simple unified theory for flow in the canopy and roughness sublayer. Boundary-Layer Meteorology, 123, 339–363.
• Bonan, G. (2019). Climate Change and Terrestrial Ecosystem Modeling. Cambridge University Press.

rsl_corrected_profile(uf_params, rsl_model, Δz_eff, ζ, z0, transport, scheme = PointValueScheme())

Return the RSL-corrected dimensionless profile F̂ = F + P, where F is the MOST profile UF.dimensionless_profile and P the roughness sublayer correction (see rsl_profile_correction, also for the arguments). All exchange coefficients, similarity scales, the bulk Richardson number, and profile recovery use this function, so that the RSL correction is applied consistently.

Since φ̂ = φ μ with μ ≤ 1, the exact corrected profile satisfies F̂ ≥ F; this bound is enforced (P ≥ 0) to guard against quadrature error, so F̂ > 0 whenever F > 0.

rsl_profile_correction(uf_params, rsl_model, Δz_eff, ζ, z0, transport, scheme = PointValueScheme())

Return the roughness sublayer correction P = F̂ - F ≥ 0 to the MOST dimensionless profile F for the given transport type, stability parameter ζ = Δz_eff / L, and discretization scheme, where F̂ is the RSL-corrected profile (see rsl_corrected_profile).

For point values, with z_c = min(Δz_eff, z_RSL) (limited to [z0, max(z_RSL, z0)]),

P = ∫_{z_c}^{z_RSL} φ(z/L) (1 - μ(z)) dz/z,

which vanishes above the RSL. For layer averages (LayerAverageScheme), P is the layer average (1/Δz_eff) ∫_{z0}^{Δz_eff} P_point(z) dz of the point-value correction, consistent with the layer-averaged MOST profile. The integrals are evaluated by Gauss-Legendre quadrature in ln z (in z for the part of the layer average that is an integral with respect to z).

Arguments

• uf_params: Universal function parameters.
• rsl_model: RSL model (LinearRSL, ExponentialRSL, or NoRoughnessSubLayer).
• Δz_eff: Effective height Δz - d above the displacement height [m].
• ζ: Stability parameter Δz_eff / L [-].
• z0: Roughness length for the transported quantity [m].
• transport: UF.MomentumTransport or UF.HeatTransport.
• scheme: PointValueScheme (default) or LayerAverageScheme.

Returns

The correction P [-], zero for NoRoughnessSubLayer and above the RSL.

NoStabilityCap()

No cap on the stability parameter: standard MOST (the default).

ConstantStabilityCap(ζ_max)

Cap on the stability parameter entering the flux-profile relations at the constant ζ_max > 0. Physick & Garratt (1995) limit z/L to 0.5 in their mesoscale model; caps between 0.5 and 2 are common in land models. The cap must be positive, so that unstable conditions are unaffected. In the flux computation, the cap is converted to the floating-point type of the parameter set. The MOST solve brackets the root for caps within its stability range ζ ≤ 100 (see surface_fluxes).

Fields

• ζ_max: Maximum stability parameter entering the flux-profile relations [-].

Examples

config = SurfaceFluxConfig(roughness, gustiness, MoistModel(), NoRoughnessSubLayer(),
    ConstantStabilityCap(0.5))

References

• Physick, W. L., & Garratt, J. R. (1995). Incorporation of a high-roughness lower boundary into a mesoscale model for studies of dry deposition over complex terrain. Boundary-Layer Meteorology, 74, 55–71.

MaxHeatFluxStabilityCap()

Cap on the stability parameter at the stability $ζ_p$ at which the MOST sensible heat flux at fixed wind speed, $H ∝ ζ / F̂_m(ζ)^3$, is maximal.

$ζ_p$ depends only on the universal functions, the solver scheme, the ratio of the effective height to the momentum roughness length (Δz_eff / z0m), and the roughness sublayer correction (if any); it has no free parameters. For point values, it satisfies $F_m(ζ_p) = 3 [φ_m(ζ_p) - φ_m(ζ_p z_{0m}/Δz_{eff})]$; for log-linear stable functions $φ_m = 1 + b ζ$, $ζ_p ≈ \ln(Δz_{eff}/z_{0m}) / (2 b)$. With the Gryanik et al. (2020) functions (point values, no RSL), $ζ_p ≈ 0.15–0.3$ for Δz_eff / z0m = 2–10 (tall canopies, forcing height a few roughness lengths above the displacement height), $≈ 0.4–0.8$ for Δz_eff / z0m = 30–1000 (grass), $≈ 1–1.2$ for Δz_eff / z0m = 3000–10⁴ (bare soil), and $≈ 1.4–1.6$ for Δz_eff / z0m = 3·10⁴–10⁵ (snow).

Up to $ζ_p$, the fluxes are those of MOST. Beyond it, the exchange coefficients are held at their values at $ζ_p$, so the heat flux increases with the surface–air temperature difference; MOST alone predicts a decreasing heat flux there, which is dynamically unstable (runaway cooling). At the cap, the local flux Richardson number $R_f = ζ_p / φ_m(ζ_p)$ is ≈ 0.09–0.22 for the Gryanik functions (increasing with Δz_eff / z0m), i.e., at or below the critical value $R_{f,cr} ≈ 0.20–0.25$ beyond which local similarity theory ceases to apply (Grachev et al. 2013); the bulk ratio $ζ_p / F_m(ζ_p)$ is ≈ 0.07–0.14 (≈ $1/(3b)$ for log-linear functions).

The cap is computed once per surface flux solve from the momentum roughness length at neutral stability, by golden-section maximization of $ζ / F̂_m(ζ)^3$ in $\log ζ$ over $ζ ∈ [10^{-2}, 20]$ with a fixed number of iterations, refined by one parabolic step (see max_heat_flux_stability). The maximization evaluates the dimensionless momentum profile 15 times per solve; when z0m depends on the friction velocity (e.g., COARE3RoughnessParams), a neutral solve for the roughness length precedes it. Since $ζ_p$ depends only on the geometry Δz_eff / z0m, the scheme, and the RSL parameters, it is constant in time for fixed roughness lengths; for such surfaces, ConstantStabilityCap(max_heat_flux_stability(...)) computed once per column gives the same fluxes, with the maximization done once.

Examples

config = SurfaceFluxConfig(roughness, gustiness, MoistModel(), NoRoughnessSubLayer(),
    MaxHeatFluxStabilityCap())

References

• Derbyshire, S. H. (1999). Boundary-layer decoupling over cold surfaces as a physical boundary-instability. Boundary-Layer Meteorology, 90, 297–325.
• van de Wiel, B. J. H., et al. (2012). The minimum wind speed for sustainable turbulence in the nocturnal boundary layer. Journal of the Atmospheric Sciences, 69, 3116–3127.
• Grachev, A. A., Andreas, E. L, Fairall, C. W., Guest, P. S., & Persson, P. O. G. (2013). The critical Richardson number and limits of applicability of local similarity theory in the stable boundary layer. Boundary-Layer Meteorology, 147, 51–82.
• Gryanik, V. M., Lüpkes, C., Grachev, A., & Sidorenko, D. (2020). New modified and extended stability functions for the stable boundary layer based on SHEBA and parametrizations of bulk transfer coefficients for climate models. Journal of the Atmospheric Sciences, 77, 2687–2716.

max_heat_flux_stability(param_set, Δz_eff, z0m, scheme = PointValueScheme(), rsl_model = NoRoughnessSubLayer())

Return the stability parameter $ζ_p$ at which $ζ / F̂_m(ζ)^3$ (the MOST sensible heat flux at fixed wind speed, up to a factor) is maximal, where F̂_m = F_m + P_m is the roughness-sublayer-corrected dimensionless momentum profile (rsl_corrected_profile).

The maximum is located by golden-section search in $\ln ζ$ over $[10^{-2}, 20]$ with a fixed number of iterations (no data-dependent control flow), followed by one parabolic (Newton) step. The step refines the value to $\approx 10^{-4}$ relative accuracy and makes $ζ_p$ a smooth function of Δz_eff, z0m, and the RSL parameters with the correct (implicit-function) derivatives under automatic differentiation. The result lies within the search interval; if the maximum lies outside it, the nearest end point is returned.

Arguments

• param_set: Parameter set.
• Δz_eff: Effective height Δz - d above the displacement height [m].
• z0m: Momentum roughness length [m].
• scheme: PointValueScheme (default) or LayerAverageScheme.
• rsl_model: Roughness sublayer model (default NoRoughnessSubLayer).

Returns

The stability parameter $ζ_p = Δz_eff / L$ of maximum heat flux [-].

See also MaxHeatFluxStabilityCap.

neutral_momentum_roughness(roughness_model, param_set, inputs, scheme)

Return the momentum roughness length z0m [m] at neutral stability. For a roughness model independent of the friction velocity (see depends_on_ustar), this is the model's roughness length; for the others (e.g., COARE3RoughnessParams), it comes from a neutral solve with compute_ustar_and_roughness.

stability_cap_value(stability_cap, param_set, inputs, scheme[, z0m])

Return the numerical value of the cap on the stability parameter (or nothing for NoStabilityCap) for the given inputs. For MaxHeatFluxStabilityCap, the momentum roughness length z0m is evaluated at neutral stability with neutral_momentum_roughness unless it is given (e.g., from a prescribed friction velocity).

with_stability_cap(inputs, param_set, scheme[, z0m])

Return inputs with the field ζ_cap set to the numerical value of the stability cap from stability_cap_value (or nothing for NoStabilityCap). The MOST solver and the prescribed-flux paths call this once per solve, so that the cap is computed once and read by capped_stability thereafter.

resolved_stability_cap(param_set, inputs, scheme)

Return the numerical value of the stability cap for inputs, or nothing: inputs.ζ_cap if set (inside the MOST solver and the prescribed-flux paths, see with_stability_cap), otherwise computed from inputs.stability_cap with stability_cap_value. NoStabilityCap and ConstantStabilityCap return their value directly; MaxHeatFluxStabilityCap runs the maximization in max_heat_flux_stability (and, when z0m depends on the friction velocity, a neutral solve for the roughness length) on each call, so callers that evaluate several quantities from the same inputs should set the cap once with with_stability_cap.

capped_stability(ζ, ζ_cap)
capped_stability(inputs, ζ)
capped_stability(param_set, inputs, scheme, ζ)

Return the stability parameter entering the flux-profile relations: min(ζ, ζ_cap), or ζ if there is no cap (ζ_cap === nothing). The second method reads the cap from inputs.ζ_cap, which is set by the MOST solver (see with_stability_cap); inputs built without a solve have ζ_cap = nothing. The third method computes the cap from inputs.stability_cap when inputs.ζ_cap is nothing (see resolved_stability_cap), so that the exchange coefficients and similarity scales computed from builder inputs agree with those of the solve.

UniversalFunctions

Universal stability and stability correction functions for the SurfaceFluxes module.

Supports the following flux-profile (ϕ, ψ, Ψ) parameterizations:

• Businger: Businger et al. (1971), Dyer (1974)
• Gryanik: Gryanik et al. (2020)
• Grachev: Grachev et al. (2007)

Both standard finite-difference (point-value) and finite-volume (layer-averaged) schemes are supported; the finite-volume scheme follows Nishizawa & Kitamura (2018).

The module also provides variance/TKE similarity functions (MomentumVariance, HeatVariance) from Panofsky et al. (1977), Wyngaard et al. (1971), and Tan et al. (2018). These are empirical surface-layer closures that are independent of the flux-profile parameterization above (see the phi(..., MomentumVariance()) / phi(..., HeatVariance()) docstrings for their range of validity).

References

• Businger, J. A., Wyngaard, J. C., Izumi, Y., & Bradley, E. F. (1971). Flux-profile relationships in the atmospheric surface layer. Journal of the Atmospheric Sciences, 28, 181–189. DOI: 10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2
• Dyer, A. J. (1974). A review of flux-profile relationships. Boundary-Layer Meteorology, 7, 363–372. DOI: 10.1007/BF00240838
• Gryanik, V. M., Lüpkes, C., Grachev, A., & Sidorenko, D. (2020). New modified and extended stability functions for the stable boundary layer based on SHEBA and parametrizations of bulk transfer coefficients for climate models. Journal of the Atmospheric Sciences, 77, 2687–2716. DOI: 10.1175/JAS-D-19-0255.1
• Grachev, A. A., Andreas, E. L., Fairall, C. W., Guest, P. S., & Persson, P. O. G. (2007). SHEBA flux–profile relationships in the stable atmospheric boundary layer. Boundary-Layer Meteorology, 124, 315–333. DOI: 10.1007/s10546-007-9177-6
• Nishizawa, S., & Kitamura, Y. (2018). A surface flux scheme based on the Monin-Obukhov similarity for finite volume models. Journal of Advances in Modeling Earth Systems, 10, 3159–3175. DOI: 10.1029/2018MS001534
• Panofsky, H. A., Tennekes, H., Lenschow, D. H., & Wyngaard, J. C. (1977). The characteristics of turbulent velocity components in the surface layer under convective conditions. Boundary-Layer Meteorology, 11, 355–361. DOI: 10.1007/BF02186086
• Wyngaard, J. C., Coté, O. R., & Izumi, Y. (1971). Local free convection, similarity, and the budgets of shear stress and heat flux. Journal of the Atmospheric Sciences, 28, 1171–1182. DOI: 10.1175/1520-0469(1971)028<1171:LFCSAT>2.0.CO;2
• Panofsky, H. A., & Dutton, J. A. (1984). Atmospheric Turbulence: Models and Methods for Engineering Applications. Wiley, New York, 397 pp.
• Tan, Z., Kaul, C. M., Pressel, K. G., Cohen, Y., Schneider, T., & Teixeira, J. (2018). An extended eddy-diffusivity mass-flux scheme for unified representation of subgrid-scale turbulence and convection. Journal of Advances in Modeling Earth Systems, 10, 770–800. DOI: 10.1002/2017MS001162

phi(p, ζ, transport)

Universal (similarity) function: the non-dimensional vertical gradient of wind shear (ϕ_m, MomentumTransport) or of temperature/scalars (ϕ_h, HeatTransport) at stability parameter ζ.

Dispatches on the parameterization type of p (BusingerParams, GryanikParams, GrachevParams) and on transport.

psi(p, ζ, transport_type)

The standard integrated stability correction function ψ(ζ). Defined as: ψ(ζ) = ∫[0 to ζ] (ϕ(0) - ϕ(x)) / x dx

This is the standard correction used in point-based Monin-Obukhov similarity theory.

Psi(p, ζ, transport_type)

The layer-averaged stability correction function Ψ(ζ). Mathematically, this is defined as: Ψ(ζ) = (1/ζ) ∫[0 to ζ] ψ(x) dx

This function is required for finite-volume models where fluxes are calculated using cell-averaged values rather than point values at the cell center.

See Nishizawa & Kitamura (2018), Eqs. 14 & 15.

BusingerParams{FT}

Parameter bundle for the Businger (1971) similarity relations. Mappings to Nishizawa & Kitamura (2018) coefficients:

• a_m, a_h: The linear coefficients for stable conditions (β in some texts).
• b_m, b_h: The coefficients γ inside the unstable sqrt/cbrt terms (e.g., (1 - γζ)).
• Pr_0: The neutral Prandtl number.

See Businger et al. (1971) and Nishizawa & Kitamura (2018).

GryanikParams{FT}

Parameter bundle for the Gryanik et al. (2020) similarity relations. These functions are designed to be valid across the entire stability range, including very stable conditions.

• a_m, b_m: Coefficients for momentum stability function (Eq. 32).
• a_h, b_h: Coefficients for heat stability function (Eq. 33).
• Pr_0: Neutral Prandtl number. (The paper recommends Pr_0 ≈ 0.98.)

Reference: Gryanik et al. (2020).

GrachevParams{FT}

Parameter bundle for the Grachev et al. (2007) similarity relations, based on SHEBA data.

• a_h, b_h, c_h: Coefficients for heat/scalar stability function (Eq. 9b). Note: c_h is the coefficient for the linear ζ term in the denominator.
• Pr_0: Neutral Prandtl number. In the original Grachev et al. (2007) derivation, this is 1.0. It is included here for structural consistency with other parameterizations.

Reference: Grachev et al. (2007).

MomentumTransport

Type selecting momentum-transfer stability functions (ϕₘ, ψₘ, Ψₘ).

HeatTransport

Type selecting heat/scalar-transfer stability functions (ϕₕ, ψₕ, Ψₕ).