Radiative Transfer

Breeze.jl integrates with RRTMGP.jl to provide radiative transfer capabilities for atmospheric simulations. The radiative transfer model computes longwave and shortwave radiative fluxes, which can be incorporated into energy tendency equations.

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 Breeze
using Breeze.AtmosphereModels
using Oceananigans.Units
using Dates
using RRTMGP

Nz = 64
λ, φ = -70.9, 42.5  # longitude, latitude
grid = RectilinearGrid(size=Nz, x=λ, y=φ, z=(0, 20kilometers),
                       topology=(Flat, Flat, Bounded))

# Thermodynamic setup
surface_temperature = 300
constants = ThermodynamicConstants()

reference_state = ReferenceState(grid, constants;
                                 surface_pressure = 101325,
                                 potential_temperature = surface_temperature)

dynamics = AnelasticDynamics(reference_state)

# Create gray radiation model
radiation = 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 position
clock = 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:

# Longwave fluxes (ZFaceFields)
ℐ_lw_up = radiation.upwelling_longwave_flux
ℐ_lw_dn = radiation.downwelling_longwave_flux

# Shortwave flux (direct beam only for non-scattering solver)
ℐ_sw = radiation.downwelling_shortwave_flux

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 = 0
solar_position = DiurnalSolarPosition(latitude = 30)

# Perpetual June solstice analog at 45°N
solar_position = DiurnalSolarPosition(latitude = 45, declination = 23.5)

# Non-Earth rotator with a 10-hour day; the simulation starts at sunrise
solar_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_solar_zenith_angle! 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.Units
using RRTMGP, NCDatasets # Required for RRTMGP lookup tables

grid = 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; the three optics flavors share the solar-position machinery.

Surface Properties

The RadiativeTransferModel model requires surface properties:

Integration with dynamics

Radiative fluxes can be used to compute heating rates for the energy equation. The radiative heating rate is computed from flux divergence:

\[F_{\mathscr{I}} = -\frac{1}{\rho cᵖᵐ} \frac{\partial \mathscr{I}_{net}}{\partial z}\]

where $\mathscr{I}_{net}$ is the net radiative flux (upwelling minus downwelling), $cᵖᵐ$ is the mixture heat capacity, and $F_{\mathscr{I}}$ is the radiative flux divergence (heating rate).

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.