Ideally, we should perform line-by-line three-dimensional (3D) RT computations to emulate as closely as possible what is happening in nature, but that is not a practical proposition. Interestingly, the practical solutions brought to bear by the almost non-intersecting spectroscopic and 3D RT communities are conceptually similar, and not only because they rely ultimately on the simpler textbook case of 1D RT, for which a plethora of computational methods exist:
To cope with spectral variability, the popular "correlated-k" methods tally across wavelengths in a given band how often a certain cross-section is encountered going up and down the absorption coefficient scale. Then one runs the 1D RT model just once for each bin in cross-section value, and uses the appropriately weighted sum of these 1D RT results to predict the outcome for the whole band.
To cope with spatial variability, one first notes that the most important and variable optical property of a 3D scene is the horizontal distribution of optical depth (OD) in each column. Therefore, tally the histogram of OD across the scene, perform a 1D RT computation for each OD bin and compute the weighted sum of these results across the bins. This is the "independent column approximation" (ICA) as applied to unresolved spatial variability in a certain domain.
The Monte Carlo ICA (McICA) approach [Barker et al., 2008] adopted in many Global Climate Models (GCMs) simply blurs the line between these spatial and spectral integrations and implements a combined random quadrature rule, primarily to deal with complications in cloud overlap across layers (e.g., "maximum/random" overlap rule). The correlated-k methodology also has difficulties when scattering is non-negligible. To be incorporated effectively into a standard multiple scattering model, one needs transmission to be expressed in the form e-<k>s where <k> is the mean extinction coefficient in the band and s = h/cosθ is the geometric path across a layer (thickness h, polar angle θ). That path can vary over at least an order-of-magnitude, thus further challenging the correlated-k technique already dealing with the extreme variability of k itself. The ultimate price is systematic underestimation of gaseous transmission.
In GCM-type radiation budget estimation, angular integration (focus on fluxes rather than radiances) is assumed from the outset, so two-stream approximations in 1D RT are widely adopted. In sharp contrast, remote sensing exploits radiance sampled usually in a small number of directions. Consequently, a full 1D RT model is required in the ICA framework and the main approximation is about neglecting the net horizontal transport inside the scene (sub-pixel 3D RT effects) and/or at its boundaries (3D RT pixel-adjacency effects).
We propose an alternative unified approach to spatial and/or spectral variability that may prove to be better suited for remote sensing applications. Davis  and Conley and Collins  independently proposed to use the same non-exponential transmission law, 1/(1+<k>s/a)a, to capture respectively the effects of spatial and spectral variability in extinction. The classic exponential law is recovered in the limit a → ∞ (no variance in the extinction coefficient k). Therefore, 1/a is a metric of relative variability in k, either spatial or spectral. The implicit assumption here is that the variability of k can be approximated by a Gamma distribution.
Davis  had scattering in mind while Conley and Collins  were thinking strictly about absorption. The latter demonstrated numerically the superiority of their parameterization over the standard exponential form for a water vapor band, ignoring scattering; in particular, the systematic bias of the exponential law is removed. The former showed that numerical (Monte Carlo) simulations of multiple scattering controlled by the non-exponential transmission law better fits ground-based oxygen A-band spectroscopy of broken and/or multi-layered cloudy skies. Moreover, Davis  reminds us of Barker et al.'s  empirical evidence that histograms of optical depth, which is proportional to the layer-averaged extinction, for a wide variety of cloud fields can be fit reasonably well by Gamma distributions.
Barker  used that evidence to design the Gamma-weighted two-stream incarnation of the ICA, which has the advantage of being analytically tractable in closed form. (It is now superseded by McICA nonetheless.) Davis  showed that the non-exponential transmission law leads to a whole new class of generalized RT models defined by RT equations in integral form that do not have obvious integro-differential counterparts.
We will survey previous research and showcase recently developed deterministic solvers for generalized 1D RT in the presence of multiple scattering based on Markov chain formalism [Xu et al., 2011]. There are both two-stream (flux-based) and N-stream (radiance-based) versions to support radiative budget and remote sensing applications respectively with highly efficient computations.
HW Barker, JAS, 53, 2289- (1996).
HW Barker et al., JAS, 53, 2304- (1996).
HW Barker et al., QJRMS, 134, 1463-(2008).
AJ Conley and WD Collins, JQSRT, 112, 1525- (2011).
AB Davis, Lecture Notes in Computational Sci. and Eng., 48, 84- (2006).
F Xu et al., Opt. Exp., 19, 946 (2011).
Supplementary URL: http://science.jpl.nasa.gov/people/ADavis/