Chapter 7: The two-stream solver

Once gas optics has produced optical depth \(\tau\), single-scattering albedo \(\omega_0\), and asymmetry parameter \(g\) for every (band, g-point, layer), a radiative-transfer solver converts them to fluxes. climt uses the Meador & Weaver (1980) two-stream formulation with \(\delta\)-Eddington scaling — the same physics as RRTMGP (Pincus et al. 2019).

The two-stream approximation

The radiative-transfer equation in a plane-parallel layer is reduced to two angular moments: the upward hemispheric flux \(F^+\) and the downward \(F^-\). For each layer this yields the diffuse reflectance \(R_\text{dif}\) and transmittance \(T_\text{dif}\). climt’s per-layer SW coefficients are computed in _sw_dif_and_source:

import inspect
from climt._components.cork.sw.kernels import _sw_dif_and_source
print(inspect.getsource(_sw_dif_and_source))

@njit
def _sw_dif_and_source(tau, w0, g, mu0):
    """Compute per-layer diffuse R/T and direct-beam source terms.

    Following RRTMGP's sw_dif_and_source (mo_rte_solver_kernels.F90).
    Uses Zdunkowski PIFM gamma coefficients.

    Args:
        tau: layer optical depth (scalar, delta-scaled)
        w0: single scattering albedo (scalar, delta-scaled)
        g: asymmetry parameter (scalar, delta-scaled)
        mu0: cosine of solar zenith angle (scalar)

    Returns:
        Rdif: diffuse reflectance
        Tdif: diffuse transmittance
        Rdir: fraction of direct beam scattered upward (source_up per unit flux)
        Tdir: fraction of direct beam scattered downward (source_dn per unit flux)
        Tnoscat: direct beam transmittance (Beer's law)
    """
    # PIFM gamma coefficients (Zdunkowski et al. 1980)
    gamma1 = (8.0 - w0 * (5.0 + 3.0 * g)) * 0.25
    gamma2 = 3.0 * (w0 * (1.0 - g)) * 0.25

    # k = sqrt(gamma1^2 - gamma2^2), floored to avoid div by 0 at conservative limit
    k = np.sqrt(max((gamma1 - gamma2) * (gamma1 + gamma2), _MIN_K))
    exp_minusktau = np.exp(-tau * k)
    exp_minus2ktau = exp_minusktau * exp_minusktau

    # RT_term: common factor (refactored for numerical stability)
    RT_term = 1.0 / (k * (1.0 + exp_minus2ktau) + gamma1 * (1.0 - exp_minus2ktau))

    # Diffuse reflectance and transmittance (M&W Eqs. 25, 26)
    Rdif = RT_term * gamma2 * (1.0 - exp_minus2ktau)
    Tdif = RT_term * 2.0 * k * exp_minusktau

    # Direct beam: unscattered transmittance
    mu0_s = max(mu0, _MIN_MU0)
    Tnoscat = np.exp(-tau / mu0_s)

    # Direct beam source terms (M&W Eqs. 14, 15)
    k_mu = k * mu0_s
    denom_dir = 1.0 - k_mu * k_mu
    if abs(denom_dir) < 1e-30:
        denom_dir = 1e-30

    RT_term_dir = w0 * RT_term / denom_dir

    gamma3 = (2.0 - 3.0 * mu0_s * g) * 0.25
    gamma4 = 1.0 - gamma3
    alpha1 = gamma1 * gamma4 + gamma2 * gamma3  # Eq. 16
    alpha2 = gamma1 * gamma3 + gamma2 * gamma4  # Eq. 17

    k_gamma3 = k * gamma3
    k_gamma4 = k * gamma4

    Rdir = RT_term_dir * (
        (1.0 - k_mu) * (alpha2 + k_gamma3)
        - (1.0 + k_mu) * (alpha2 - k_gamma3) * exp_minus2ktau
        - 2.0 * (k_gamma3 - alpha2 * k_mu) * exp_minusktau * Tnoscat
    )

    Tdir = -RT_term_dir * (
        (1.0 + k_mu) * (alpha1 + k_gamma4) * Tnoscat
        - (1.0 - k_mu) * (alpha1 - k_gamma4) * exp_minus2ktau * Tnoscat
        - 2.0 * (k_gamma4 + alpha1 * k_mu) * exp_minusktau
    )

    # Energy clamping (Robin Hogan / Peter Ukkonen)
    Rdir = max(0.0, min(Rdir, 1.0 - Tnoscat))
    Tdir = max(0.0, min(Tdir, 1.0 - Tnoscat - Rdir))

    return Rdif, Tdif, Rdir, Tdir, Tnoscat

Figure 1: Figure 7.1 — Diffuse reflectance \(R\) and transmittance \(T\) vs. layer optical depth for conservative scattering (\(\omega_0 = 1\), \(g = 0.5\)). A thin layer mostly transmits; a thick layer mostly reflects.

δ-Eddington scaling

Real cloud and aerosol phase functions have a strong forward peak that the two-stream cannot resolve. The \(\delta\)-Eddington transformation absorbs the forward peak into an effective reduction of optical depth:

\[\tau' = \tau(1 - \omega_0 g^2), \quad \omega'_0 = \frac{\omega_0(1-g^2)}{1-\omega_0 g^2}, \quad g' = \frac{g}{1+g}.\]

climt applies this inside _delta_scale before every two-stream solve.

The adding method

With per-layer \((R, T)\) known, the adding method combines layers bottom-up into a column albedo and flux profile. The _adding function in climt/_components/cork/sw/kernels.py propagates both diffuse and direct-beam components in one downward-then-upward sweep.

Longwave transport

The LW solver is simpler — no scattering, just absorption and emission. A diffusivity factor \(D = 1.66\) converts vertical \(\tau\) to effective-diffuse \(\tau\). The compiled kernel loops over all (band, g-point, column) triples:

import inspect
from climt._components.cork.lw.kernels import _lw_transport_kernel
print(inspect.getsource(_lw_transport_kernel))

@njit(parallel=True)
def _lw_transport_kernel(
    tau, planck_source, surface_source, emissivity, weights,
    up_band, down_band, up_broad, down_broad,
    diag_trans, diag_up_gpt, diag_dn_gpt, want_diag, diffusivity_factor,
):
    """Consolidated multi-band, multi-g-point LW transport.

    Loops over columns in parallel; for each (band, g-point) runs the up/down
    diffusivity sweeps and accumulates weighted fluxes into up_band/down_band
    inside the compiled kernel. Accumulation order (g ascending, then b
    ascending for broadband) matches the original python loops bit-for-bit.
    """
    nband, ngpt, nlev, ncol = tau.shape
    for i in prange(ncol):
        for k in range(nlev + 1):
            up_broad[k, i] = 0.0
            down_broad[k, i] = 0.0
        for b in range(nband):
            for k in range(nlev + 1):
                up_band[b, k, i] = 0.0
                down_band[b, k, i] = 0.0
            for g in range(ngpt):
                w = weights[b, g]
                # Upward sweep: surface -> TOA
                up_prev = emissivity[b, i] * surface_source[b, g, i]
                up_band[b, 0, i] += w * up_prev
                if want_diag != 0:
                    diag_up_gpt[b, g, 0, i] = w * up_prev
                for k in range(nlev):
                    trans = np.exp(-diffusivity_factor * tau[b, g, k, i])
                    up_cur = up_prev * trans + planck_source[b, g, k, i] * (1.0 - trans)
                    up_band[b, k + 1, i] += w * up_cur
                    if want_diag != 0:
                        diag_trans[b, g, k, i] = trans
                        diag_up_gpt[b, g, k + 1, i] = w * up_cur
                    up_prev = up_cur
                # Downward sweep: TOA -> surface (dn_prev starts at 0 = TOA BC)
                dn_prev = 0.0
                if want_diag != 0:
                    diag_dn_gpt[b, g, nlev, i] = 0.0
                for k in range(nlev - 1, -1, -1):
                    trans = np.exp(-diffusivity_factor * tau[b, g, k, i])
                    dn_cur = dn_prev * trans + planck_source[b, g, k, i] * (1.0 - trans)
                    down_band[b, k, i] += w * dn_cur
                    if want_diag != 0:
                        diag_dn_gpt[b, g, k, i] = w * dn_cur
                    dn_prev = dn_cur
            for k in range(nlev + 1):
                up_broad[k, i] += up_band[b, k, i]
                down_broad[k, i] += down_band[b, k, i]

The diffusivity factor \(D = 1.66\) is the Elsasser approximation; it converts the exact hemispheric integral \(\int_0^1 e^{-\tau/\mu}\,d\mu\) to \(e^{-D\tau}\).

Further reading

• Meador and Weaver (1980) — the two-stream coefficients and their derivation.
• Pincus et al. (2019) — RRTMGP, the modern successor; discusses the adding method and \(\delta\)-Eddington in production use.

References

Meador, W. E., and W. R. Weaver. 1980. “Two-Stream Approximations to Radiative Transfer in Planetary Atmospheres: A Unified Description and Results.” Journal of the Atmospheric Sciences 37 (3): 630–43. https://doi.org/10.1175/1520-0469(1980)037<0630:TSATRT>2.0.CO;2.

Pincus, R., E. J. Mlawer, and J. S. Delamere. 2019. “Balancing Accuracy, Efficiency, and Flexibility in Radiation Calculations for Dynamical Models.” Journal of Advances in Modeling Earth Systems 11: 3074–89. https://doi.org/10.1029/2019MS001621.