Radiative Transfer

Breeze.jl computes longwave and shortwave radiative fluxes with a RadiativeTransferModel, whose flux divergence enters the energy tendency of an AtmosphereModel. Two radiation backends are available as package extensions, selected by the optics passed to the constructor:

• RRTMGP.jl solves gray-atmosphere (GrayOptics), clear-sky (ClearSkyOptics), and all-sky (AllSkyOptics) radiation with its correlated-$k$ lookup tables; load it with using RRTMGP, ClimaComms, NCDatasets.
• NumericalRadiation.jl solves clear-sky and all-sky radiation with the ecCKD gas optics (EcCKDOptics) in one kernel per column, on CPUs and GPUs; load it with using NumericalRadiation: NumericalRadiation and using NCDatasets (NumericalRadiation exports its own ThermodynamicConstants, so it is loaded qualified next to Breeze). See ecCKD radiation with NumericalRadiation.jl below.

Both backends fill the same four flux fields and the same flux divergence, so a simulation switches between them by changing the optics alone.

Gray Atmosphere Radiation

The simplest radiative transfer option is gray atmosphere radiation, which uses the optical thickness parameterization from Schneider (2004) and O'Gorman and Schneider (2008). This approximation treats the atmosphere as having a single effective absorption coefficient rather than computing full spectral radiation.

Basic Usage

To use gray radiation in a Breeze simulation, create a RadiativeTransferModel model with the GrayOptics optics flavor and pass it to the AtmosphereModel constructor:

using Breezeusing Breeze.AtmosphereModelsusing Oceananigans.Unitsusing Datesusing RRTMGPNz = 64λ, φ = -70.9, 42.5  # longitude, latitudegrid = RectilinearGrid(size=Nz, x=λ, y=φ, z=(0, 20kilometers),                       topology=(Flat, Flat, Bounded))# Thermodynamic setupsurface_temperature = 300constants = ThermodynamicConstants()reference_state = ReferenceState(grid, constants;                                 base_pressure = 101325,                                 potential_temperature = surface_temperature)dynamics = AnelasticDynamics(reference_state)# Create gray radiation modelradiation = RadiativeTransferModel(grid, GrayOptics(), constants;                                   surface_temperature,                                   surface_emissivity = 0.98,                                   surface_albedo = 0.1,                                   solar_constant = 1361) # W/m²# Create atmosphere model with DateTime clock for solar positionclock = Clock(time=DateTime(2024, 9, 27, 16, 0, 0))model = AtmosphereModel(grid; clock, dynamics, radiation)

AtmosphereModel{CPU, RectilinearGrid}(time = 2024-09-27T16:00:00, iteration = 0)
├── grid: 1×1×64 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── dynamics: AnelasticDynamics(p₀=101325.0, θ₀=300.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float64}
├── timestepper: SSPRungeKutta3
├── advection scheme: 
│   ├── momentum: Centered(order=2)
│   ├── ρθ: Centered(order=2)
│   └── ρqᵛ: Centered(order=2)
├── forcing: @NamedTuple{ρu::Returns{Float64}, ρv::Returns{Float64}, ρw::Returns{Float64}, ρθ::Returns{Float64}, ρqᵛ::Returns{Float64}, ρE::Returns{Float64}}
├── tracers: ()
├── coriolis: Nothing
└── microphysics: Nothing

When a DateTime clock is used, the cosine of the solar zenith angle is computed automatically from the time and grid location (longitude and latitude). See Solar zenith angle below for the full set of options.

Gray Radiation Model

The RadiativeTransferModel model computes:

• Longwave radiation: Both upwelling and downwelling thermal radiation using RRTMGP's two-stream solver
• Shortwave radiation: Direct beam solar radiation

The gray atmosphere optical thickness for longwave follows the parameterization by O'Gorman and Schneider (2008),

\[τ_{lw} = α \frac{Δp}{p_0} \left[ f_l + 4 (1 - f_l) \left(\frac{p}{p_0}\right)^3 \right] \left[ τ_e + (τ_p - τ_e) \sin^2 φ \right]\]

where $φ$ is latitude and $α$, $f_l$, $τ_e$, and $τ_p$ are empirical parameters.

For shortwave:

\[τ_{sw} = 2 τ_0 \frac{Δp}{p_0} \frac{p}{p_0}\]

where $τ_0 = 0.22$ is the shortwave optical depth parameter.

The above two expressions are identical to those in the RRTMGP documentation.

Radiative Fluxes

After running set!, the radiative fluxes are available from the radiation model as four ZFaceFields, one per stream. All three optics flavors expose the same four fields, and all of them follow the sign convention positive upward, so downwelling fluxes are stored as negative numbers:

# Longwaveℐ_lw_up = radiation.upwelling_longwave_flux     # ≥ 0ℐ_lw_dn = radiation.downwelling_longwave_flux   # ≤ 0# Shortwaveℐ_sw_up = radiation.upwelling_shortwave_flux    # ≥ 0ℐ_sw_dn = radiation.downwelling_shortwave_flux  # ≤ 0

For clear-sky and all-sky optics, downwelling_shortwave_flux is the total downwelling shortwave (direct beam plus diffuse), and upwelling_shortwave_flux carries the sunlight scattered back upward by air and clouds together with the fraction reflected by the surface.

Solar zenith angle

The cosine of the solar zenith angle $\cos(θ_z)$ controls the magnitude of the top-of-atmosphere shortwave flux and the slant-path optical depth through the column. Breeze provides one keyword — solar_position — that selects how $\cos(θ_z)$ is determined on each radiation update. It takes any subtype of AbstractSolarPosition. The two concrete subtypes cover the common cases:

Time-varying apparent sun

ApparentSolarPosition accepts two optional keyword arguments:

• coordinate: an explicit (longitude, latitude) tuple in degrees, or nothing (the default) to read $(λ, φ)$ from the grid's coordinates. For single-column grids the grid's (x, y) is interpreted as (λ, φ).
• epoch: a DateTime anchor for floating-point model clocks. The model time in seconds is added to epoch to obtain the absolute DateTime. With a DateTime clock, epoch is ignored.

# Today's default: DateTime clock, λ/φ inferred from the grid.solar_position = ApparentSolarPosition()# Float clock + epoch — useful on lat–lon / curvilinear grids where clock# precision matters but you want full per-column zenith.solar_position = ApparentSolarPosition(epoch = DateTime(2024, 1, 1))# Pin a 3D simulation to a specific observer (overrides per-column lat/lon).solar_position = ApparentSolarPosition(coordinate = (-70.9, 42.5),                                       epoch      = DateTime(2024, 1, 1))

Idealized diurnal cycle

For idealized studies where you want a clean 24-hour cycle (or any rotation period) without orbital mechanics, calendar dates, or equation-of-time corrections, use DiurnalSolarPosition:

# Perpetual equinox at 30°N — 24-hour day, noon at t = 0solar_position = DiurnalSolarPosition(latitude = 30)# Perpetual June solstice analog at 45°Nsolar_position = DiurnalSolarPosition(latitude = 45, declination = 23.5)# Non-Earth rotator with a 10-hour day; the simulation starts at sunrisesolar_position = DiurnalSolarPosition(latitude = 0,                                      day_length  = 10 * 3600,                                      noon_offset = 5  * 3600)

The cosine of the solar zenith angle is computed analytically each radiation update from the model clock:

\[\cos(θ_z) = \sin(φ) \sin(δ) + \cos(φ) \cos(δ) \cos(ω), \qquad ω = \frac{2π}{T_d} (t - t_{\text{noon}}),\]

where $φ$ is the (fixed) latitude, $δ$ is the (fixed) declination, $T_d$ is the day length, and $t_{\text{noon}}$ is the simulation time at which local noon occurs. There is no annual cycle, no equation of time, and no longitude — by design.

Required: numeric clock

DiurnalSolarPosition reads model.clock.time as "seconds since the start of the simulation" and computes the hour angle from it. A DateTime clock is rejected with an actionable error, because pairing one with this idealized cycle would be semantically ambiguous (the calendar would be ignored). Use Clock(time = 0.0) (or any numeric Clock) when constructing the AtmosphereModel.

Choosing noon_offset

The default noon_offset = 0 places local noon at t = 0 — the sun is overhead at the start of the simulation. To start at sunrise instead, set noon_offset = day_length / 4; to start at midnight, set noon_offset = day_length / 2. The cosine function is naturally periodic, so no mod is required and the cycle simply repeats every day_length seconds.

Gray-optics latitude

Unlike FixedCosineZenith — where latitude for the gray τ comes from the grid — DiurnalSolarPosition carries its own latitude and uses it for both the diurnal cos(θ_z) calculation and the gray-optics latitude-dependent τ. This keeps the two consistent in idealized single-column setups where the grid's y might be a placeholder.

Constant cos(θ_z)

For idealized studies where you want a single, time-independent solar geometry — typical in RCE intercomparisons — use FixedCosineZenith:

radiation = RadiativeTransferModel(grid, GrayOptics(), constants;                                   surface_temperature  = 300,                                   surface_albedo       = 0.1,                                   solar_constant       = 1361,    # W/m²                                   solar_position       = FixedCosineZenith(0.5))# A floating-point clock works fine — no epoch is required.clock = Clock(time = 0.0)model = AtmosphereModel(grid; clock, dynamics, radiation)

The cosine of the zenith angle is written into the RRTMGP boundary-condition array once at construction and never recomputed; the per-step update_cos_zenith! call becomes a no-op.

Choosing a value

Common choices for $\cos(θ_z)$ in idealized work:

Interaction with solar_constant

The top-of-atmosphere downward shortwave flux is solar_constant * cos_zenith, so solar_constant and FixedCosineZenith together control both:

• the TOA SW magnitude (their product), and
• the slant-path absorption (which depends on $\cos(θ_z)$ alone, through $\exp(-τ/\cos(θ_z))$).

Note that scaling solar_constant and scaling cos_zenith are not equivalent for the shortwave heating profile. The TOA flux changes the same way, but the shape of the absorption with height does not. If your study cares about the vertical structure of SW heating (it usually does), pick $\cos(θ_z)$ to match the slant path you actually want, then choose solar_constant to set the magnitude.

For example, to model diurnal-mean conditions with the full solar constant spread over the day's path:

solar_position = FixedCosineZenith(0.5)solar_constant = 1361 / 2      # because TOA SW = solar_constant * cos_zenith

Latitude for gray optics

The gray-optics longitudinal-mean optical thickness depends on latitude through $τ_e$ and $τ_p$. With FixedCosineZenith you specify only $\cos(θ_z)$ — Breeze still reads the latitude needed for the gray τ from the grid's coordinates (or from coordinate, in 3D setups where you pass an explicit position via ApparentSolarPosition). This means it's fine to combine a fixed zenith with a single-column grid located at a particular latitude: τ_e/τ_p will reflect that latitude, while the zenith stays pinned.

Clear-sky Full-spectrum Radiation

For more accurate radiative transfer calculations, use the ClearSkyOptics optics flavor which computes full-spectrum gas optics using RRTMGP's lookup tables:

using Breeze, Oceananigans.Unitsusing RRTMGP, NCDatasets # Required for RRTMGP lookup tablesgrid = RectilinearGrid(; size=16, x=0, y=45, z=(0, 10kilometers),                       topology=(Flat, Flat, Bounded))constants = ThermodynamicConstants()radiation = RadiativeTransferModel(grid, ClearSkyOptics(), constants;                                   surface_temperature = 300,                                   surface_emissivity = 0.98,                                   surface_albedo = 0.1,                                   background_atmosphere = BackgroundAtmosphere(CO₂ = 400e-6))

RadiativeTransferModel
├── solar_constant: 1361.0 W m⁻²
├── solar_position: ApparentSolarPosition(coordinate=(0.0, 45.0), epoch=<from clock>)
├── surface_temperature: ConstantField(300.0) K
├── surface_emissivity: ConstantField(0.98)
├── direct_surface_albedo: ConstantField(0.1)
└── diffuse_surface_albedo: ConstantField(0.1)

The BackgroundAtmosphere struct specifies volume mixing ratios for radiatively active gases (CO₂, CH₄, N₂O, O₃, etc.). Water vapor is computed from the model's prognostic moisture field.

Clear-sky and all-sky models accept the same solar_position keyword as gray-optics; all optics flavors share the solar-position machinery.

ecCKD radiation with NumericalRadiation.jl

EcCKDOptics selects the ecCKD correlated-$k$ gas optics of Hogan and Matricardi (2022), solved by NumericalRadiation.jl. The gas optics tables ship with the ecRad data that NumericalRadiation downloads as an artifact, and are read from netCDF files, so both NumericalRadiation and NCDatasets must be loaded. NumericalRadiation exports its own ThermodynamicConstants; loading it qualified avoids the clash with Breeze's:

using Breeze, Oceananigans.Unitsusing NCDatasetsusing NumericalRadiation: NumericalRadiationgrid = RectilinearGrid(; size=16, x=0, y=45, z=(0, 10kilometers),                       topology=(Flat, Flat, Bounded))constants = ThermodynamicConstants()radiation = RadiativeTransferModel(grid, EcCKDOptics(), constants;                                   surface_temperature = 300,                                   surface_emissivity = 0.98,                                   surface_albedo = 0.1,                                   background_atmosphere = BackgroundAtmosphere(CO₂ = 400e-6),                                   overlying_atmosphere = OverlyingAtmosphereColumn(top = 65kilometers, layers = 40))

RadiativeTransferModel
├── solar_constant: 1361.0 W m⁻²
├── solar_position: ApparentSolarPosition(coordinate=(0.0, 45.0), epoch=<from clock>)
├── surface_temperature: ConstantField(300.0) K
├── surface_emissivity: ConstantField(0.98)
├── direct_surface_albedo: ConstantField(0.1)
├── liquid_effective_radius: ConstantRadiusParticles{Float64}(1.0e-5)
├── ice_effective_radius: ConstantRadiusParticles{Float64}(3.0e-5)
├── diffuse_surface_albedo: ConstantField(0.1)
├── optics: EcCKDOptics with 32 longwave and 32 shortwave g-points, clear sky
└── overlying_atmosphere: 40 layers from 10000.0 m to 65000.0 m

The keyword arguments are those of the RRTMGP clear-sky and all-sky constructors (surface properties, solar_position, solar_constant, schedule, background_atmosphere, effective radii) plus overlying_atmosphere. The default EcCKDOptics() uses the :climate_32x32 reference model with 32 longwave and 32 shortwave g-points; EcCKDOptics(:climate_64x64) doubles the spectral resolution, and a (longwave = path, shortwave = path) pair selects ecCKD definition files of your own. On a GPU() grid the tables and the solve live on the device.

Overlying atmosphere

A limited-area grid rarely reaches the top of the atmosphere. Rather than treating the grid top as the top of the atmosphere, an ecCKD model continues the radiation column above it through a OverlyingAtmosphereColumn: by default 40 geometrically stretched layers from the grid top to 65 km, whose temperature follows the U.S. Standard Atmosphere (standard_atmosphere_temperature) anchored to the temperature on the grid's top face and relaxing to the standard profile over a 1 km blending_height. Pressure is hydrostatic above the grid, the specific humidity follows an idealized profile floored at the stratospheric value, and the ozone mole fraction follows BackgroundAtmosphere.O₃ — a number or a function of height, evaluated above the grid as well — unless ozone_mole_fraction overrides it. Extension layers are clear.

The extension is what gives the downwelling longwave flux on the grid's top face its stratospheric contribution (a few W m⁻² for a 20 km grid) and lets the stratospheric ozone absorb ultraviolet sunlight before it reaches the grid, instead of depositing that absorption in the grid's top cells. Pass overlying_atmosphere = nothing to stop the column at the grid top, which is the right choice for comparisons with the RRTMGP models (whose column ends there) or for a grid that already reaches the mesosphere.

# Stop at the grid topradiation = RadiativeTransferModel(grid, EcCKDOptics(), constants;                                   surface_temperature = 300, surface_albedo = 0.1,                                   overlying_atmosphere = nothing)# A taller, finer extension with a constant stratospheric ozone mole fractionextension = OverlyingAtmosphereColumn(top = 80kilometers, layers = 60, ozone_mole_fraction = 2e-6)

Gases

The ecCKD reference models tabulate eight gases. Two vary in space: water vapor, taken from the model's prognostic moisture, and ozone, from BackgroundAtmosphere.O₃. Five are well-mixed and read from the BackgroundAtmosphere: CO₂, CH₄, N₂O, CFC₁₁ and CFC₁₂. Nitrogen and oxygen form the dry air, whose absorption the tables fold in. Gas amounts follow the "dry" column convention the ecCKD tables were derived with: the dry-air amount of a layer is its total mass over the dry molar mass, ρ Δz / mᵈ (Δp / (g mᵈ) for a hydrostatic layer), the water vapor ρ qᵛ Δz / mᵛ is counted on top of it, and the well-mixed gases are their mole fractions times the dry-air amount. The molar masses mᵈ and mᵛ and the gravitational acceleration g are those of the ThermodynamicConstants passed to RadiativeTransferModel; the Stefan–Boltzmann constant of the gray Planck source (used by g-points without a Planck source table) is carried by NumericalRadiation's gas optics model. The remaining gases of BackgroundAtmosphere (CO, NO₂, CFC₂₂, CCl₄, CF₄ and the HFCs) are not tabulated and must be left at zero; a nonzero value throws an ArgumentError rather than being silently ignored.

Clouds

All-sky ecCKD radiation adds the cloud liquid and ice diagnosed by the model's microphysics to the layer optics through single-scattering tables mapped onto the same g-points as the gas optics, selected with CloudScatteringTables:

optics = EcCKDOptics(clouds = CloudScatteringTables(liquid = :mie_droplet,                                                    ice = :baum_general_habit_mixture))radiation = RadiativeTransferModel(grid, optics, constants;                                   surface_temperature = 300, surface_albedo = 0.1,                                   liquid_effective_radius = ConstantRadiusParticles(10e-6),                                   ice_effective_radius = ConstantRadiusParticles(30e-6))

The two selectors name the Mie droplet and Baum general-habit-mixture ice tables shipped with the ecRad data; a path to a netCDF file in the same format works in their place. The tables are evaluated once, at construction, at the constant effective radii of liquid_effective_radius and ice_effective_radius (only ConstantRadiusParticles for now). The shortwave folds the cloud scattering into the two-stream solution with delta-Eddington scaling; the longwave adds the cloud absorption only and neglects longwave cloud scattering (ecRad's do_lw_cloud_scattering = false setting, as RRTMGP's no-scattering longwave solver does). With clouds = nothing (the default) the sky is clear whatever the microphysics holds.

Sign and units conventions

NumericalRadiation solves each column top-down with longwave and shortwave fluxes positive in their own direction of travel. Breeze converts on the way back to the grid, so an ecCKD model exposes the same four ZFaceFields as the RRTMGP models, in W m⁻², all positive upward with downwelling fluxes stored as negative numbers, and the same flux_divergence in W m⁻³ (see Integration with dynamics). On the grid's top face the downwelling shortwave equals -solar_constant * cos_zenith only without a column extension; with one it is reduced by the absorption above the grid. The downwelling longwave on the top face is zero without an extension and the emission of the atmosphere above with one.

Coexistence with RRTMGP

The two extensions load side by side (using RRTMGP, ClimaComms, NCDatasets together with using NumericalRadiation: NumericalRadiation), and GrayOptics, ClearSkyOptics and AllSkyOptics remain RRTMGP models. Because both backends fill the same fields, one script can build an RRTMGP model and an ecCKD model on two AtmosphereModels with identical initial conditions and difference their fluxes; the single column radiation example does exactly that. On the same 20 km column, with overlying_atmosphere = nothing so that both stop at the grid top, the two agree to 2 W m⁻² in the clear-sky boundary fluxes, to 0.1 W m⁻² in the 2×CO₂ forcing and to 0.2 K day⁻¹ (root mean square) in the clear-sky heating rates, the spread expected between two independent correlated-$k$ models; they part at cloud top, where the cloud optics differ in kind.

Surface Properties

The RadiativeTransferModel model requires surface properties:

Integration with dynamics

Radiative fluxes provide a volumetric energy source from flux convergence,

\[Q_{\mathscr{I}} = -\frac{\partial \mathscr{I}_{net}}{\partial z},\]

where $Q_{\mathscr{I}}$ has units W/m³. For AnelasticDynamics, holding moisture phase composition fixed, the corresponding temperature tendency is

\[\left. \frac{\partial T}{\partial t} \right|_{rad} = \frac{Q_{\mathscr{I}}}{ρᵣ cᵖᵐ}.\]

Here $ρᵣ$ is the reference density, $cᵖᵐ$ is the mixture heat capacity, and $\mathscr{I}_{net}$ is the net upward radiative flux. Because Breeze stores downwelling components as negative values, the four signed streams are added:

\[\mathscr{I}_{net} = \mathscr{I}_{lw}^{↑} + \mathscr{I}_{lw}^{↓} + \mathscr{I}_{sw}^{↑} + \mathscr{I}_{sw}^{↓} .\]

Despite its historical name, radiation.flux_divergence stores $Q_{\mathscr{I}}$: minus the vertical derivative of the signed flux sum, evaluated at cell centers from the two bounding cell faces. The static-energy formulation adds this W/m³ source directly to its energy-density tendency. The liquid-ice potential-temperature formulation adds $Q_{\mathscr{I}} / (\hat{c}ᵖ \hat{Π})$ to its potential-temperature-density tendency, with the vapor-equivalent heat capacity $\hat{c}ᵖ$ and Exner function $\hat{Π}$ (see Liquid-ice potential temperature). Omitting the upwelling shortwave would bias either tendency wherever shortwave radiation scatters or reflects.

Architecture Support

The radiative transfer implementation supports both CPU and GPU architectures. The column-based RRTMGP solver is called from Oceananigans' field data arrays with appropriate data layout conversions. The ecCKD backend runs on the grid's architecture without leaving the device. It can solve the columns in batches of whole grid rows (column_batches, as for RRTMGP; one solve of the whole grid by default): per batch, a layer pass reads the model's fields directly and stores each layer's gas-optics stencil, gas amounts, cloud water paths and interface temperatures, one work item per layer, and a column kernel streams the g points through them and accumulates the fluxes directly into the radiation model's flux fields. Only one batch of per-column storage is allocated, so setting column_batches caps the radiation's memory on large grids, at some cost in speed on a GPU.