Spherical harmonic radiance code

Overview

The spherical harmonic radiance code solves the monochromatic radiative transfer equation for radiances rather than integrated fluxes, expanding the angular dependence of the radiation field in spherical harmonics. This provides greater angular resolution than the two-stream approach, at higher computational cost, and is particularly suited to applications requiring accurate radiance fields (e.g. remote sensing simulations).

The monochromatic equation of transfer is:

where \(I\) is the radiance, \(k^{(s)}\) and \(k^{(a)}\) are the scattering and absorption coefficients, \(P(n', n)\) is the phase function, and \(j\) is the emission source term.

The phase function is expanded in Legendre polynomials:

and may be \(\delta\)-rescaled using the standard prescription \(k^{(s)} \to (1-f)k^{(s)}\) before further treatment.

Fundamentals of solving for radiances

Assuming a plane-parallel atmosphere, the radiance field is expanded in spherical harmonics \(Y_l^m(n)\):

where the direct solar beam \(I_\odot\) is kept separate. Using the orthogonality of the \(Y_l^m\) and introducing optical depth \(\tau\) and single-scattering albedo \(\omega = k^{(s)} / (k^{(s)} + k^{(a)})\), the equation of transfer reduces to a coupled system:

where \(s_l = 1 - \omega g_l\) and the Clebsch–Gordan coefficients are:

The atmosphere is divided into \(N\) homogeneous layers with optical thicknesses \(\tau_i\), \(i = 1, \ldots, N\). The full solution in each layer is the sum of a complementary function and a particular integral.

The complementary function

The complementary function consists of exponentials \(H_{lm}(\mu) e^{\tau/\mu}\) for \(\mu \in \mathbb{R}\). For fixed azimuthal order \(m\), the permissible values of \(\mu\) are defined by an eigenvalue problem. Truncating the spherical harmonic expansion at odd order \(L\), the recurrence relations:

where \(L' = L\) if \(m\) is even and \(L' = L+1\) if \(m\) is odd. Defining \(K_{lm} = \sqrt{s_l} H_{lm}\), this becomes a symmetric tridiagonal eigenvalue problem:

with non-zero entries \(C_{l-1,l} = c^+_{l-1,m}/\sqrt{s_l s_{l-1}}\).

Numerical stability: For \(|\mu| > 1\) the recurrence has growing solutions triggered by rounding errors. The recurrence is therefore always applied in the direction of decreasing \(l\). When scattering is nearly conservative, intermediate terms can overflow; the recurrence is recast in \(H'_{lm} = \sigma^{-l} H_{lm}\) where \(\sigma = 1/\max(1, 2\mu^{-1})\) to prevent this.

Halving the eigenproblem: Since eigenvalues come in \(\pm\mu_k\) pairs, the even and odd components of the eigenvector decouple, halving the size of the matrix to be diagonalised. The eigenvalues of the reduced matrix have diagonal elements \(d_j = E^2_{2j-1} + E^2_{2j}\) and sub-diagonal elements \(e_j = E_{2j-2} E_{2j-1}\).

The complementary function in the \(i\)-th layer is written using only negative exponentials to avoid overflow:

Conservative scattering (\(\omega \to 1\)) causes the matrix \(C\) to become singular for \(m = 0\). This is handled by artificially reducing \(\omega\) slightly to avoid ill-conditioning.

Particular integrals

Thermal radiation

In the infra-red, the equation of transfer is reformulated in terms of differential radiances \(I' = I - B\). For a Planckian source varying linearly across a layer:

with \(I_{i,lm} = 0\) otherwise, where \(\Delta B_i\) is the difference in the Planck function across layer \(i\).

For quadratic variation of the Planckian:

Small optical depths cause ill-conditioning. This is resolved by adding to the particular integral a solution of the homogeneous system that exhibits the same singularity as \(\tau \to 0\):

where \(V_{kl} = K_{klm}/\sqrt{s_l}\) and the coefficients \(u^\pm_k = \mp q_0 s_1 V_{k1}\) are determined by matching the singularity structure.

Solar particular integral

Using the convention \(\mu_0 = -\cos\theta_\odot\), the direct solar beam in layer \(i\) is:

A particular integral of the form \(I_{ilm}(\tau) = Z_{ilm} e^{-\tau/\mu_0}\) is sought, giving the recurrence:

with solution:

where \(\gamma = I_\odot(\Delta_{i-1}) Y_{L'+1}^{m*}(n_\odot) / V_{i,L'+1,m}(\mu_0)\) enforces the truncation condition, and \(V(\mu_0)\) satisfies the associated homogeneous recurrence.

Ill-conditioning arises when \(\mu_0\) approaches an eigenvalue \(\mu_k\). This is removed by subtracting a weighted multiple of the corresponding eigensolution, applied smoothly via a weighting that depends on the separation \(|\mu_0 - \mu_k|\) rather than an if-test.

Boundary conditions

Interior boundaries

Continuity of the radiance field at each interior level requires:

Upper boundary

At the top of the atmosphere, the downward radiance is specified as \(I(n) = I^{(0)}(n)\) for \(n \in \Omega^-\). Since the full boundary condition cannot be imposed in a truncated system, Marshak's conditions ¹ are used — the inner product of the residual with odd-parity harmonics is set to zero:

where:

For \(l + l'\) even, \(\kappa_{ll'm} = \frac{1}{2} \delta_{ll'}\). For \(l + l'\) odd, \(\kappa_{ll'm}\) is evaluated analytically using associated Legendre polynomial identities:

Surface boundary conditions

The surface is characterised by a bidirectional reflectance distribution function (BRDF) \(\gamma_r(n, n')\), so that:

For a Lambertian surface, \(\gamma_r = \alpha/\pi\) where \(\alpha\) is the albedo.

The BRDF is expanded in a double spherical harmonic series:

Imposing rotational symmetry, reflectional symmetry, reciprocity (\(\gamma_r(n,n') = \gamma_r(n',n)\)), and the reality of \(\gamma_r\) reduces the coefficients to real, symmetric quantities \(\Psi_{ll'm}\) with \(\Psi_{ll'-m} = \Psi_{ll'm}\).

The full surface boundary condition including thermal emission is:

where \(B^*\) is the isotropic Planckian radiance at the surface temperature, arising from Kirchhoff's law.

After applying Marshak's procedure and expanding, the boundary condition becomes:

where \(\Xi\), \(\Phi\), and \(\Lambda\) are precomputed geometric coupling matrices that depend only on the BRDF expansion coefficients \(\rho_j\) and the \(\kappa_{ll'm}\) integrals.

The BRDF of the ocean surface

The ocean surface BRDF is not available from a single direct reference ². The radiance code itself is used to compute the radiance within the ocean, with special upper boundary conditions to handle refraction. Key optical properties:

Rayleigh scattering in water:

where \(\lambda_0 = 550\,\text{nm}\), \(K_R = 0.93\) (pure water) or \(1.21\) (sea water).

Particulate scattering ³ is related to chlorophyll concentration \(C\) (mg m\(^{-3}\)):

Absorption ^{4 5}:

Refraction at the surface follows Snell's law (\(\sin\theta_i = n\sin\theta_t\), \(n \approx 1.34\)), with the Fresnel reflection coefficient for unpolarised light:

Total internal reflection occurs for upward rays beyond the critical angle (\(r_{wa} = 1\)). The boundary condition in the ocean is:

Numerical implementation

Since \(I \in \mathbb{R}\), coefficients satisfy \(I_{l,-m} = (-1)^m I_{lm}^*\), so only \(m \geq 0\) need be stored. The complex \(e^{im\phi}\) dependence factors out as \(I_{lm} = C_{lm} e^{-im\phi_\odot}\), and the full radiance field reduces to:

reducing storage requirements by roughly a factor of two compared to a naive complex implementation.

Increasing the speed of computation

The source function (TMS) technique

Computing radiances by spherical harmonics alone converges slowly because high orders are needed to represent singly scattered radiation. The TMS method ⁶ separates single and multiple scattering: the full (unrescaled) phase function is used for single scattering, while the rescaled truncated phase function handles multiple scattering.

The diffuse radiance equation under rescaling becomes:

The algorithm proceeds as follows:

1. Solve the rescaled, truncated equation using spherical harmonics to obtain \(\hat{I}_T\) and \(\hat{I}_\odot\).
2. Integrate along a ray using \(\hat{I}_T\) as the source function for multiply scattered radiation.
3. Compute the unrescaled single-scattering contribution from the true solar beam \(I_\odot\) separately.

The final integration along a ray between optical depths \(\hat{\Delta}^-\) and \(\hat{\Delta}^+\) gives:

where \(A_{ilm}\) and \(B_i\) are layer contributions evaluated by integrating exponential terms analytically.

Ill-conditioning near eigenvalue crossings

When \(\mu \to \tilde{\mu}\) (an eigenvalue or \(\mu_0\)), geometric factors of the form \(\tilde{\mu}/(\tilde{\mu} - \mu)\) become ill-conditioned. This is regularised using a perturbation \(\eta\):

where \(\varepsilon\) is machine epsilon, so that:

This introduces errors of \(O(\sqrt{\varepsilon})\) only in a neighbourhood of the singularity.

Fast solution of the linear equations

A more efficient algorithm for the core linear system — not yet fully implemented — reduces the dominant operation count from \(O(18N^3 L)\) (banded solver with partial pivoting) to \(O(6N^3 L)\).

Block structure

For fixed azimuthal order \(m\), retaining \(2N\) polar orders, the amplitudes at the top and bottom of each layer are assembled into a block matrix system. Collecting alternate (even/odd) eigenvector components as \(W_{sk}\) and \(U_{sk}\) respectively, the orthogonality relations give:

Marshak's condition at the top boundary becomes a block equation involving the matrix \(M_{sp} = \tilde{M}_{2s+1-m,\, 2p-m}\).

Forward recurrence

After block Gaussian elimination, the system reduces to a recurrence of the form:

where the matrices \(A_n\), \(B_n\), \(C_n\), \(D_n\) are related to the eigenvector blocks and layer transmission factors \(\theta_k = e^{-\bar{\tau}/\mu_k}\).

Key simplifications arise from the relations \(R = P^{-1}\), \(S = Q^{-1}\), and:

which reduce the cost of computing \(P\) and \(R\) from \(O(4N^3)\) to \(O(4N^2)\) multiplications.

Back substitution

After forward elimination, back substitution proceeds as:

Only matrices \(A\), \(B\), \(X\) need be retained at each level, reducing memory by a factor of three relative to the banded solver.