Add `BoundaryLayerAtmosphere` — a bespoke prognostic boundary-layer atmosphere model

[glwagner]

Motivation

PrescribedAtmosphere (JRA55 / ECCO / ERA5 / OSPapa) and full prognostic atmospheres (Breeze, SpeedyWeather) cover the two ends of the cost/fidelity spectrum, but nothing in between. The CheapAML formulation (Deremble, Wienders & Dewar, MWR 2013, doi:10.1175/MWR-D-11-00254.1) fills that gap: a 1.5-D single-layer marine boundary-layer model that prognoses near-surface air temperature $T$ and specific humidity $q$ only, taking winds and downwelling radiation from a reanalysis. Much cheaper than a full atmosphere, but avoids the SST-decoupling pathology of fully-prescribed forcing in long ocean integrations.

Scope of this PR

• In scope: a new prognostic atmosphere model (struct, time-stepping, advection, source terms, optional diagnostic precipitation) that plugs into the existing ESM coupling slot.
• Not in scope: any change to surface flux computation. $\mathcal{Q}^T$, $\mathcal{Q}^v$, $J^v$, momentum stress, and the assembled net surface heat flux $\Sigma \mathcal{Q}^{ao}$ are all produced by NumericalEarth's existing similarity-theory and ocean-flux-assembly machinery (SimilarityTheoryFluxes, compute_atmosphere_ocean_fluxes!, _assemble_net_ocean_fluxes!, COARE3-equivalent). BLA reads those fluxes as inputs to its tracer tendencies; energy conservation across the air–sea interface is automatic because the heat flux at the BL bottom is $\Sigma \mathcal{Q}^{ao}$.

Proposed design: BoundaryLayerAtmosphere <: AbstractModel

The new type is a bespoke model class (following the ClimaSeaIce.SeaIceModel precedent), not a wrapper around PrescribedAtmosphere. It mirrors PrescribedAtmosphere's outward field names so the two are interchangeable from the ESM's point of view; only tracers differs in kind — prescribed FieldTimeSeries for PrescribedAtmosphere, prognostic Fields for BoundaryLayerAtmosphere. Default physics are CheapAML; the type's name describes what the model is, leaving room for future variants under the same name.

It implements the standard Oceananigans AbstractModel interface — time_step!(model::BoundaryLayerAtmosphere, Δt; callbacks) and update_state!(model) — so the existing ESM time_step!(coupled_model, Δt) loop dispatches to it directly with no new indirection.

struct BoundaryLayerAtmosphere{...} <: AbstractModel{TS, Arch}
    architecture
    grid                                 # 2D, topology (·, ·, Flat)
    clock
    # Prognostic state
    tracers                              # (T, q) — Fields at (C, C, Nothing)
    # Prescribed forcing inputs (same names as PrescribedAtmosphere)
    velocities                           # PrescribedVelocityFields backed by FieldTimeSeries
    pressure                             # FieldTimeSeries
    downwelling_radiation                # TwoBandDownwellingRadiation (ℐꜜˢʷ, ℐꜜˡʷ)
    # Precipitation regime — flat layout, selector decides what runs:
    #   precipitation === nothing                 → rain/snow are FieldTimeSeries (prescribed)
    #   precipitation isa DiagnosticPrecipitation → rain/snow are writable Fields,
    #                                               populated each step by the LSP+CP kernel
    rain                                 # FieldTimeSeries | Field
    snow                                 # FieldTimeSeries | Field
    precipitation                        # Nothing | DiagnosticPrecipitation
    auxiliary_freshwater_flux            # rivers/icebergs, or nothing
    thermodynamics_parameters
    surface_layer_height
    boundary_layer_height                # h_bℓ
    # Local buffer of AO interface fluxes — populated by
    # update_net_fluxes!(coupled_model, ::BoundaryLayerAtmosphere) at the end of each
    # outer update_state!. The next time_step! reads from this buffer (one-step lagged).
    # Same pattern Breeze uses (`ext/NumericalEarthBreezeExt/breeze_atmosphere_interface.jl:91-140`).
    surface_fluxes                       # NamedTuple of CenterFields
    # Numerics — plug-and-play from Oceananigans
    advection                            # default WENO()
    closure                              # default nothing (WENO carries enough dissipation)
    timestepper                          # default :RungeKutta3
    forcing                              # tracer source terms
    # Physics knobs
    longwave_emissivity, moisture_entrainment_factor, air_density
end

A small parameter struct selects diagnostic precipitation and carries its tunables:

struct DiagnosticPrecipitation{FT}
    large_scale_timescale       :: FT   # τ₁, default 40 h
    convective_timescale        :: FT   # τ₂, default 6 h
    reference_vertical_velocity :: FT   # w₀, default 7.5e-6 m/s
    convective_threshold        :: FT   # q* threshold, default 2e-4 kg/kg
    freezing_temperature        :: FT   # rain/snow split, default 273.15 K
end

Two constructors: a primary one taking forcing inputs explicitly, and a one-line convenience that lifts everything off an existing PrescribedAtmosphere:

prescribed = JRA55PrescribedAtmosphere(arch)

# Diagnostic precipitation (default)
atmos = BoundaryLayerAtmosphere(prescribed; grid = ocean_grid)

# Prescribed precipitation pass-through (use the dataset's rain/snow)
atmos = BoundaryLayerAtmosphere(prescribed; grid = ocean_grid,
                                rain          = prescribed.freshwater_flux.rain,
                                snow          = prescribed.freshwater_flux.snow,
                                precipitation = nothing)

model = AtmosphereOceanModel(atmos, ocean)

Users can also mix sources (JRA55 winds + ERA5 radiation + OSPapa pressure) by calling the primary constructor directly.

API note. This flattens the precipitation pair out of the nested freshwater_flux NamedTuple that PrescribedAtmosphere currently uses. Adopting the same flat layout in PrescribedAtmosphere (along with a precipitation = nothing marker) would keep the two types structurally identical. Worth a follow-up issue.

Coupling sequence (one ESM step)

time_step!(coupled_model, Δt):
    time_step!(sea_ice, Δt)
    time_step!(ocean, Δt)
    time_step!(atmosphere, Δt)              # standard AbstractModel dispatch — Oceananigans RK3:
                                            #   compute_tendencies! reads model.surface_fluxes
                                            #   (last update's AO fluxes), prescribed radiation,
                                            #   and assembles BL-top + BL-bottom + entrainment
                                            #   + relaxation source terms
                                            #   ... time-stepper integrates T, q
                                            #   ... LSP+CP kernel populates rain/snow if
                                            #       precipitation isa DiagnosticPrecipitation
    tick!(coupled_model.clock, Δt)
    update_state!(coupled_model):
        interpolate_state! for each component
        compute_atmosphere_ocean_fluxes!    # uses BLA's T, q via the exchanger
        update_net_fluxes!(model, atmos::BoundaryLayerAtmosphere)
                                            # writes ao.sensible_heat, ao.latent_heat,
                                            # ao.water_vapor and the radiative components
                                            # of ΣQao into model.surface_fluxes
        update_net_fluxes! for ocean, sea_ice

Forward-Euler-on-couplings with one-step-lagged surface fluxes inside BLA's RHS — the same lagging the existing slab/ocean components use (src/Oceans/slab_ocean.jl:140-148).

Why bespoke AbstractModel (not HFSM, not a wrapper)

• It's a model: prognostic state, tendencies, time-stepper.
• Mirroring PrescribedAtmosphere's API (same field names) means downstream consumers don't care which kind they hold.
• Plugs into Oceananigans' machinery: any AbstractAdvectionScheme, any AbstractTurbulenceClosure, any AbstractTimeStepper, plus Simulation/output-writers/checkpointer/CFL-wizard for free. Time-stepping reuses Oceananigans.TimeSteppers.TimeStepper(:RungeKutta3, …) — no custom integrator, and no new ESM-side dispatch. :SplitRungeKutta3 is not applicable since the model has no implicit substep.
• Source terms (BL-bottom and BL-top heat fluxes, BL-top moisture entrainment, Newtonian relaxation) become DiscreteForcings on the tracers; closures capture references to model.surface_fluxes and the prescribed downwelling radiation FTS.

Turbulence

No turbulence closure required for this model. With WENO advection, no explicit Laplacian diffusion is needed either; default closure = nothing. Users wanting to experiment with HorizontalScalarDiffusivity etc. just swap the closure keyword.

Formulation summary (CheapAML, paper Eqs. 1, 4–6, 9–13)

Notation follows the NumericalEarth notation appendix: heat fluxes use $\mathcal{Q}$, vapor mass fluxes use $J^v$, radiative intensities use $\mathscr{I}$ with $\downarrow / \uparrow$ direction modifiers and $\mathrm{sw}/\mathrm{lw}$ band labels.

Prognostic equations for the two BL tracers:

$$\frac{\partial T}{\partial t} + \nabla \cdot (\mathbf{u} T) = \frac{\Sigma \mathcal{Q}^{ao} - \mathcal{Q}^{\mathrm{top}}}{\rho^{at} c^{pd} h_{b\ell}} + \nabla \cdot (K_h \nabla T) - \lambda (T - T_c)$$

$$\frac{\partial q}{\partial t} + \nabla \cdot (\mathbf{u} q) = \frac{J^v - J^{v,\mathrm{top}}}{\rho^{at} h_{b\ell}} + \nabla \cdot (K_h \nabla q) - \lambda (q - q_c)$$

BL-bottom heat flux ($\Sigma \mathcal{Q}^{ao}$): by conservation across the air–sea interface, the upward heat flux at the BL base equals the net surface heat flux NumericalEarth assembles for the ocean (_assemble_net_ocean_fluxes!:125):

$$\Sigma \mathcal{Q}^{ao} = (\mathscr{I}_\uparrow^{lw} + \mathcal{Q}^T + \mathcal{Q}^v)(1 - \aleph) + \mathscr{I}_a^{lw} + \mathcal{Q}^{ss}$$

BLA reads $\Sigma \mathcal{Q}^{ao}$ from its surface_fluxes buffer; it does not recompute components.

BL-top heat flux ($\mathcal{Q}^{\mathrm{top}}$): unique to BLA — radiative balance at the PBL top:

$$\mathcal{Q}^{\mathrm{top}} = -\mathscr{I}_\downarrow^{sw} + \frac{1}{2} \epsilon \sigma T^4 + \mathcal{Q}^v$$

where $-\mathscr{I}_\downarrow^{sw}$ accounts for shortwave passing through the BL, $\frac{1}{2}\epsilon \sigma T^4$ is the BL's upward longwave emission, and $\mathcal{Q}^v$ is the latent heat carried up by water vapor that condenses at the PBL top.

BL-bottom moisture flux ($J^v$): the same turbulent water-vapor mass flux NumericalEarth produces at atmosphere_ocean_fluxes.jl:163. Read directly from the surface_fluxes buffer.

BL-top moisture flux ($J^{v,\mathrm{top}}$): moisture entrainment across the PBL top:

$$J^{v,\mathrm{top}} = \mu \rho^{at} C^v |\mathbf{u}| q, \quad \mu = 0.25$$

The atmospheric tracer $T$ is the BL temperature, $\rho^{at}$ is air density, $c^{pd}$ is dry-air heat capacity, $\epsilon$ is BL longwave emissivity, $\sigma$ is the Stefan–Boltzmann constant, and $C^v$ is the moisture-transfer coefficient (paper's $C_{de}$).

Diagnostic precipitation (paper Eqs. 12–13; runs only when precipitation isa DiagnosticPrecipitation; populates rain/snow, partitioned by the freezing temperature):

$$w = -h_{b\ell} (\partial_x u + \partial_y v)$$

$$\mathrm{LSP} = \max \left( \rho^{at} h_{b\ell} \frac{q - 0.7 q^{sat}}{\tau_1} \left[\frac{\max(w, 0)}{w_0}\right]^2, 0 \right)$$

$$\mathrm{CP} = \max \left( \rho^{at} h_{b\ell} \frac{q - 0.9 q^{sat}}{\tau_2}, 0 \right) \quad \text{for } q > q^*_{\mathrm{thresh}}$$

with parameters from the DiagnosticPrecipitation struct.

Implementation checklist

• New src/Atmospheres/BoundaryLayerAtmosphere/ submodule:
  • boundary_layer_atmosphere.jl — struct, DiagnosticPrecipitation parameter struct, primary + PrescribedAtmosphere convenience constructors, set!, prognostic_state/restore_prognostic_state!, the AbstractModel contract methods (architecture, prognostic_fields, tracernames, etc.).
  • boundary_layer_atmosphere_interface.jl — ESM glue: thermodynamics_parameters/surface_layer_height/boundary_layer_height accessors, ComponentExchanger, interpolate_state!, update_net_fluxes!(coupled_model, ::BoundaryLayerAtmosphere) populating model.surface_fluxes.
  • boundary_layer_atmosphere_time_step.jl — time_step!(::BoundaryLayerAtmosphere, Δt; callbacks) (standard AbstractModel dispatch — no new ESM-side function), compute_tendencies!, update_state! (advances FTS for prescribed precipitation, calls the LSP+CP kernel for diagnostic).
  • boundary_layer_atmosphere_forcings.jl — DiscreteForcing constructors. The bottom-flux closure reads $\Sigma \mathcal{Q}^{ao}$ and $J^v$ from model.surface_fluxes; the top-flux closure computes BL-top radiative balance; entrainment and relaxation closures are local.
• Wire into src/Atmospheres/Atmospheres.jl; re-export BoundaryLayerAtmosphere and DiagnosticPrecipitation from src/NumericalEarth.jl.
• test/test_boundary_layer_atmosphere.jl — synthetic-FTS unit test + JRA55 integration test (both diagnostic and prescribed-precipitation paths); checkpoint round-trip; GPU.

Risk-ordered sequencing: (1) skeleton model implementing the AbstractModel contract with zero tendencies (verify ESM accepts it and its time_step! is dispatched), (2) WENO advection by prescribed velocity, (3) bottom-flux forcing reading $\Sigma \mathcal{Q}^{ao}$, (4) top-flux forcing, (5) entrainment, (6) relaxation, (7) prescribed-precipitation pass-through, (8) diagnostic precipitation, (9) checkpointing, (10) GPU.

Open design questions

1. Land relaxation reference $s_c$: not derivable automatically — user passes a relaxation = (T_c, q_c, λ_T, λ_q) NamedTuple, typically with $T_c$, $q_c$ being the prescribed-source FieldTimeSeries. Default relaxation = nothing (no relaxation).
2. Velocity staggering: WENO needs $u$ at (Face, Center, Nothing) and $v$ at (Center, Face, Nothing), but PrescribedAtmosphere.velocities is at (Center, Center, Nothing). The model allocates face-located fields and re-stages from prescribed centers in update_state!. Two extra 2D fields' worth of memory.
3. Cross-grid model: v1 requires model.grid === exchange_grid. The forcing FieldTimeSeries may live on a different grid (interpolated each step). v2 lifts the restriction (relates to Allow for pairwise exchange grids #147).
4. surface_fluxes buffer composition: which fluxes BLA mirrors. At minimum it needs $\mathcal{Q}^T$, $\mathcal{Q}^v$, $J^v$ for BL-bottom heat/moisture and the components NumericalEarth needs to assemble $\Sigma \mathcal{Q}^{ao}$ ($\mathscr{I}_\uparrow^{lw}$, $\mathscr{I}_a^{lw}$, $\mathcal{Q}^{ss}$, $\aleph$). Could alternatively mirror the assembled $\Sigma \mathcal{Q}^{ao}$ from net_fluxes.ocean.T (multiplying back by $\rho^{oc} c^{oc}$) — fewer fields, but couples BLA to a temperature-flux convention.
5. Flat rain/snow vs nested freshwater_flux: BLA exposes top-level rain, snow, precipitation. PrescribedAtmosphere currently nests rain/snow inside freshwater_flux. Opening a follow-up to flatten PrescribedAtmosphere to the same shape would make the two types structurally identical.

References

• Deremble, B., Wienders, N., & Dewar, W. K. (2013). CheapAML: A Simple, Atmospheric Boundary Layer Model for Use in Ocean-Only Model Calculations. Mon. Wea. Rev., 141, 809–821. doi:10.1175/MWR-D-11-00254.1
• Oceananigans model interface documentation
• MITgcm reference implementation: MITgcm/pkg/cheapaml/