In this section, we summarize the three two-stream models and the benchmark DISORT model with 16 streams. These algorithms are well documented in papers by Toon et al. (1989), Briegleb and Light (2007), Jin and Stamnes (1994), and Stamnes et al. (1988). Readers interested in detailed mathematical derivations should refer to those papers. We only include their key equations to illustrate the difference among two-stream models for discussion purposes.

In this work, we focus on the performance of two-stream algorithms for pure snow simulations. The inputs for these three models are the same: single-scattering properties (SSPs, i.e., single-scattering albedo ϖ, asymmetry factor g, extinction coefficient σ_ext) of snow determined by snow grain radius r, snow depth, solar zenith angle θ, solar incident flux, and the albedo of underlying ground (assuming Lambertian reflectance of 0.25 for all wavelengths). A delta transform is applied to fundamental input optical variables for all simulations (Eqs. 2a–2c, Wiscombe and Warren, 1980).

In snow, photon scattering occurs at the air–ice interface, and the absorption of photons occurs within the ice crystal. The most important factor that determines snow shortwave properties is the ratio of total surface area to total mass of snow grains, also known as “the specific surface area” (e.g., Matzl and Schneebeli, 2006, 2010). The specific surface area (β) can be converted to a radiatively effective snow grain radius r:

where ρ_ice is the density of pure ice, 917 kg m⁻³. Assuming the grains are spherical, the SSPs of snow can thus be computed using Mie theory (Wiscombe, 1980) and ice optical constants (Warren and Brandt, 2008). In nature, snow grains are not spherical, and many studies have been carried out to quantify the accuracy of such spherical representations (Grenfell and Warren, 1999; Neshyba et al., 2003; Grenfell et al., 2005). In recent years, more research has been done to evaluate the impact of grain shape on snow shortwave properties (Dang et al., 2016; He et al., 2017, 2018a, b), and they show that nonspherical snow grain shapes mainly alter the asymmetry factor. Dang et al. (2016) also point out that the solar properties of a snowpack consisting of nonspherical ice grains can be mimicked by a snowpack consisting of spherical grains with a smaller grain size by factors up to 2.4. In this work, we still assume the snow grains are spherical, and this assumption does not qualitatively alter our evaluation of the radiative transfer algorithms.

The input SSPs of snow grains are computed using Mie theory at a fine spectral resolution for a wide range of ice effective radius r from 10 to 3000 µm that covers the possible range of grain radius for snow on Earth (Flanner et al., 2007). The same spectral SSPs were also used to derive the band-averaged SSPs of snow used in SNICAR. Note Briegleb and Light (2007) refer to SSPs as inherent optical properties.

In climate modeling, snow albedo computation at a fine spectral resolution is expensive and unnecessary. Instead of computing spectrally resolved snow albedo, wider-band solar properties are more practical. For example, CESM and E3SM aggregate the narrow RRTMG bands used for the atmospheric radiative transfer simulation into visible (0.2–0.7 µm) and near-IR (0.7–5 µm) bands. The land model and sea ice model thus receive visible and near-IR fluxes as the upper boundary condition, and return the corresponding visible and near-IR albedos to the atmosphere model. In practice, these bands are also partitioned into direct and diffuse components. Therefore, a practical two-stream algorithm should be able to simulate the direct visible, diffuse visible, direct near-IR, and diffuse near-IR albedos and absorptions of snow accurately.

The band albedo α is an irradiance-weighted average of the spectral albedo α(λ):

In this work, we use the spectral irradiance F(λ) generated by the atmospheric DISORT-based Shortwave Narrowband Model (SWNB2) (Zender et al., 1997; Zender, 1999) for typical clear-sky and cloudy-sky conditions of midlatitude winter as shown in Fig. 1a. The total clear-sky down-welling surface flux at different solar zenith angles are also given in Fig. 1b.

Figure 1Spectral and total down-welling solar flux at surface computed using SWNB2 for (a) standard clear-sky and cloudy-sky atmospheric profiles of midlatitude winter assuming solar zenith angle is 60^∘ at the top of the atmosphere, and for (b) standard clear-sky profiles of midlatitude and sub-Arctic winter with different incident solar zenith angles.

Download

It has been pointed out in previous studies that the two-stream approximations become poor as solar zenith angle approaches 90^∘ (e.g., Wiscombe, 1977; Warren, 1982). As shown in Figs. 3 and 4, all three two-stream models underestimate the direct snow albedo for large solar zenith angles. In the visible band, when the snow grain size is small, the error in direct albedo is almost negligible (Fig. 3); while as snow ages and snow grains become larger, the error increases yet remains low if the snow is deep (Fig. 4). In the near-IR range, the biases of albedo are also larger for larger snow grain radii. For a given snow size, the magnitudes of such biases are almost independent of snow depth and mainly determined by the solar zenith angle. In general, the errors of all-wave direct albedo are mostly contributed by the errors of near-IR albedo, especially for optically thick snowpacks (i.e., semi-infinite), because the errors of direct albedo in the visible range are negligible compared with those in the near-IR range. To improve the performance of two-stream algorithms, we develop a parameterization that corrects the underestimated near-IR snow albedo at large zenith angles.

Figure 11 shows the direct near-IR albedo and fractional absorption of 2 m thick snowpacks consisting of grains with radii of 100 and 1000 µm, computed using two-stream algorithms and 16-stream DISORT. For solar zenith angles >75^∘, two-stream models underestimate snow albedo and overestimate solar absorption within the snowpack, mostly in the top 2 cm of snow, and the differences among the three two-stream models are small. In Sect. 5, we have shown that dEdd–AD produces the most accurate snow albedo in general. With anticipated wide application of dEdd–AD, we develop the following parameterization to adjust its low biases in computed near-IR direct albedo.

Figure 11(a) Direct near-IR snow albedo and (b) near-IR fractional absorption by top 2 cm snow of a 2 m thick snowpack, for solar zenith angles larger than 70^∘ and snow grain radii of 100 and 1000 µm. (c) The ratios of near-IR albedo computed using dEdd–AD compared to those computed using 16-stream DISORT for different solar zenith angles. These ratios are parameterized as linear functions of the logarithmic of snow grain radius. The slopes and y intercepts are shown in (d). The black dashed curves in (c, d) are fitting values computed using parameterization discussed in Sect. 5.

Download

We define and compute $R_{{\mathrm{75}}_{+}}$ as the ratio of direct semi-infinite near-IR albedo computed using 16-stream DISORT (α_16-DISORT) to that computed using dEdd–AD (α_dEdd-AD), for solar zenith angle >75^∘. This ratio is shown in Fig. 11c and can be parameterized as a function of snow grain radius (r, in meters) and the cosine of incident solar zenith angle (μ₀), as shown in Fig. 11c:

where coefficients c₁ and c₀ are polynomial functions of μ₀, as shown in Fig. 11d:

Since two-stream models always underestimate snow albedo, $R_{{\mathrm{75}}_{+}}$ always exceeds 1 (Fig. 11c). We can then adjust the direct near-IR snow albedo (α_dEdd-AD) and direct near-IR solar absorption (Fabs_dEdd-AD) by snow computed using dEdd–AD with ratio $R_{{\mathrm{75}}_{+}}$:

where F_nir is the direct near-IR flux. This adjustment reduces the error of near-IR albedo from negative 2 %–10 % to within ±0.5 % for solar zenith angles larger than 75^∘, and for grain radii ranging from 30 to 1500 µm (Fig. 12). Errors in broadband direct albedo are therefore also reduced to <0.01. The direct near-IR flux absorbed by the snowpack decreases after applying this adjustment.

Figure 12Error in semi-infinite snow albedo computed using dEdd–AD before (a, b, c) and after (d, e, f) incorporating corrections for direct near-IR albedo, for different solar zenith angles and snow grain radii.

Download

When the solar zenith angle exceeds 75^∘, our model adjusts the computed direct near-IR albedo α_dEdd−AD by the ratio $R_{{\mathrm{75}}_{+}}$ following Eqs. (12)–(14a) and reduces direct near-IR absorption following Eq. (14b). If snow is divided into multiple layers, we assume all decreased near-IR absorption (second term on the right-hand side, Eq. 14b) is confined within the top layer. This assumption is fairly accurate for the near-IR band since most absorption occurs at the surface of the snowpack (Figs. 10 and 11). As discussed previously, this parameterization is developed based on albedo computed using dEdd–AD. For models that do not use dEdd–AD but SNICAR and 2SD, the same adjustment still applies given the small differences of near-IR direct albedo computed using two-stream models (Fig. 11). For models that adopt other radiative transfer algorithms it is best for the developers to examine their model against a benchmark model such as 16-stream DISORT or two-stream models discussed in this work before applying this correction.

Although the errors of direct near-IR albedos are large for large solar zenith angles, the absolute error in reflected shortwave flux is small (Figs. 7 and 8) as the down-welling solar flux reaches snowpack and decreases as solar zenith angle increases (Fig. 1b). However, such small biases in flux can be important for high latitudes where the solar zenith angle is large for many days in late winter and early spring.

ESMs often use band-averaged SSPs of snow and aerosols for computational efficiency, rather than using brute-force integration of spectral solar properties across each band (per Eq. 11). In addition to using different radiative transfer approximations, SNICAR and dEdd–AD also adopt different methods to derive the band-averaged SSPs of snow for different band schemes.

In SNICAR, snow solar properties are computed for five bands: one visible band (0.3–0.7 µm) and four near-IR bands (0.7–1, 1–1.2, 1.2–1.5, and 1.5–5 µm). The solar properties of four subdivided near-IR bands are combined by fixed ratios to compute the direct and diffuse near-IR snow properties. These two sets of ratios are derived offline based on the incident solar spectra typical of midlatitude winter for clear- and cloudy-sky conditions (Fig. 1a).

The band-averaged SSPs of snow grains are computed following the Chandrasekhar mean approach (Thomas and Stamnes, 1999, their Eq. 9.27; Flanner et al., 2007). Specifically, spectral SSPs of snow grains are weighted into bands according to surface incident solar flux typical of midlatitude winter for clear- and cloudy-sky conditions. In addition, the single-scattering albedo ϖ(λ) of ice grains is also weighted by the hemispheric albedo α(λ) of an optically thick snowpack:

Two sets of snow band-averaged SSPs are generated for all grain radii, suitable for direct and diffuse light. For each modeling step and band, SNICAR is called twice to compute the direct and diffuse snow solar properties.

In dEdd–AD, the snow-covered sea ice properties are computed for three bands: one visible band (0.3–07 µm) and two near-IR bands (0.7–1.19 and 1.19–5 µm). The solar proprieties of these two near-IR bands are combined using ratios w_nir1 and w_nir2 for 0.7–1.19 and 1.19–5 µm, depending on the fraction of direct near-IR flux f_nidr:

The band SSPs of snow are derived by integrating the spectral SSPs and the spectral surface solar irradiance measured in the Arctic under mostly clear sky.

In addition, the band-averaged single-scattering albedo $\mathit{\varpi }\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$ is also increased to $\mathit{\varpi }{\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)}^{\prime }$ until the band albedo computed using averaged SSPs matches the band albedo $\stackrel{\mathrm{‾}}{\mathit{\alpha }}$ within 0.0001, where $\stackrel{\mathrm{‾}}{\mathit{\alpha }}$ is

dEdd–AD adopts this single set of band SSPs for both direct and diffuse computations. In practice, the physical snow grain radius r is adjusted to a radiatively equivalent radius r_eqv based on the fraction of direct flux in the near-IR band (f_nidr):

This r_eqv and the corresponding snow SSPs are then used in the radiative transfer calculation. The computed direct and diffuse solar properties alone are less accurate, while the combined all-sky broadband solar properties agree with SNICAR (Briegleb and Light, 2007). As a result, for each modeling step and band, the dEdd–AD radiative transfer subroutine is called only once to compute both the direct and diffuse snow solar properties simultaneously.

SNICAR and dEdd–AD also use different approaches to avoid numerical singularities. In SNICAR, singularities occur when the denominator of term $C_n^{\pm }$ in Eq. (3) equals zero (i.e., ${\mathit{\gamma }}^{\mathrm{2}}-\mathrm{1}/{\mathit{\mu }}_{\mathrm{0}}^{\mathrm{2}}=\mathrm{0}$), where γ is determined by the approximation method and SSPs of snow, and μ₀ is the cosine of the solar zenith angle (Eqs. 23 and 24, Toon et al., 1989). When such a singularity is detected, SNICAR will shift μ₀ by +0.02 or −0.02 to obtain physically realistic radiative properties. In the dEdd–AD algorithm, singularities arise only when μ₀=0 (Eq. 4). Therefore, in practice, for μ₀<0.01, dEdd–AD computes the sea ice solar properties for μ₀=0.01 to avoid unphysical results.

Based on the intercomparison of three two-stream algorithms and their implementations in ESMs, we formulated the following surface shortwave radiative transfer recommendations for an accurate, fast, and consistent treatment for snow on land, land ice, and sea ice in ESMs.

First, the two-stream delta-Eddington adding–doubling algorithm by Briegleb and Light (2007) is unsurpassed as a radiative transfer core. The evaluation in Sect. 5 shows that this algorithm produces the least error for snow albedo and solar absorption within snowpack, especially under overcast skies. This algorithm applies well to both uniformly refractive media such as snow on land, and to nonuniformly refractive media, such as bare, snow-covered, and ponded sea ice and bare and snow-covered land ice. Numerical singularities occur only rarely (when μ₀=0) and are easily avoided in model implementations. Among the three two-stream algorithms discussed here, dEdd–AD is also the most efficient one as it takes only two-thirds of the time of SNICAR and 2SD to compute solar properties of multilayer snowpacks.

Second, any two-stream cryospheric radiative transfer model can incorporate the parameterization described in Sect. 6 to adjust the low bias of direct near-IR snow albedo and high bias of direct near-IR solar absorption in snow, for solar zenith angles larger than 75^∘. These biases are persistent across all two-stream algorithms discussed in this work, and should be corrected for snow-covered surfaces. Alternatively, adopting a four-stream approximation would reduce or eliminate such biases, though at considerable expense in computational efficiency.

Third, in a cryospheric radiative transfer model, one should prefer physically based parameterizations that are extensible and convergent (e.g., with increasing spectral resolution) for the band-averaged SSPs and size distribution of snow. Although the treatments used in SNICAR and dEdd–AD are both practical since they both reproduce the narrowband solar properties with carefully derived band-averaged inputs as discussed in Sect. 7, the snow treatment used in SNICAR is more physically based and reproducible since it does not rely on subjective adjustment and empirical coefficients as used in dEdd–AD. Specifically, the empirical adjustment to snow grain radius implemented in dEdd–AD may not always produce compensating errors. For example, in snow containing light-absorbing impurities such adjustment may also lead to biases in aerosol absorption since the albedo reduction caused by light-absorbing particles does not linearly depend on snow grain radius (Dang et al., 2015). For further model development incorporating nonspherical snow grain shapes (Dang et al., 2016; He et al., 2018a, b), such adjustment on grain radius may fail as well. Moreover, SNICAR computes the snow properties for four near-IR bands, which helps capture the spectral variation in albedo (Fig. 2) and therefore better represents near-IR solar properties. It is also worth noting that unlike the radiative core of dEdd–AD, SNICAR is actively maintained, with numerous modifications and updates in the past decade (e.g., Flanner et al., 2012; He et al., 2018b). Snow radiative treatments that follow SNICAR conventions for SSPs may take advantage of these updates. Note that any radiative core that follows SNICAR SSP conventions must be called twice to compute diffuse and direct solar properties.

Fourth, a surface cryospheric radiative transfer model should flexibly accommodate coupled simulations with distinct atmospheric and surface spectral grids. Both the five-band scheme used in SNICAR and the three-band scheme used in dEdd–AD separate the visible from near-IR spectrum at 0.7 µm. This boundary aligns with the Community Atmospheric Model's original radiation bands (CAM; Neale et al., 2010), though not with the widely used Rapid Radiative Transfer Model (RRTMG; Iacono et al., 2008), which places 0.7 µm squarely in the middle of a spectral band. A mismatch in spectral boundaries between atmospheric and surface radiative transfer schemes can require an ESM to unphysically apportion energy from the straddled spectral bin when coupling fluxes between surface and atmosphere. The spectral grids of surface and atmosphere radiation need not be identical so long as the coarser grid shares spectral boundaries with the finer grid. In practice maintaining a portable cryospheric radiative module such as SNICAR requires a complex offline toolchain (Mie solver, spectral refractive indices for air, water, ice, and aerosols, spectral solar insolation for clear and cloudy skies) to compute, integrate, and rebin SSPs. Aligned spectral boundaries between surface and atmosphere would simplify the development of efficient and accurate radiative transfer for the coupled Earth system.

Last, it is important to note that, although we only examine the performance of the dEdd–AD for pure snow in this work, this algorithm can be applied to the surface solar calculation of all cryospheric components with or without light-absorbing particles present. First, Briegleb and Light (2007) proved its accuracy for simulating ponded and bare sea ice solar properties against observations and a Monte Carlo radiation model. Second, In CESM and E3SM, the radiative transfer simulation of snow on land ice is carried out by SNICAR with prescribed land ice albedo. Adopting the dEdd–AD radiative core in SNICAR will permit these ESMs to couple the snow and land ice as a nonuniformly refractive column for more accurate solar computations since bare, snow-covered, and ponded land ice is physically similar to bare, snow-covered, and ponded sea ice, and the latter is already treated well by the dEdd–AD radiative transfer core. Third, adding light-absorbing particles in snow will not change our results qualitatively. Both dEdd–AD and SNICAR simulate the impact of light-absorbing particles (black carbon and dust) on snow and/or sea ice using self-consistent particle SSPs that follow the SNICAR convention (e.g., Flanner et al., 2007; Holland et al., 2012). These particles are assumed to be either internally or externally mixed with snow crystals; the combined SSPs of mixtures (e.g., Appendix A of Dang et al., 2015) are then used as the inputs for radiative transfer calculation. The adoption of the dEdd–AD radiative transfer algorithm in SNICAR, and the implementation of SNICAR snow SSPs in dEdd–AD enables a consistent simulation of the radiative effects of light-absorbing particles in the cryosphere across ESM components.

In summary, this intercomparison and evaluation has shown multiple ways that the solar properties of cryospheric surfaces can be improved in the current generation of ESMs. We have merged these findings into a hybrid model SNICAR-AD, which is primarily composed of the radiative transfer scheme of dEdd–AD, five-band snow–aerosol SSPs of SNICAR, and the parameterization to correct for snow albedo biases when solar zenith angle exceeds 75^∘. This hybrid model can be applied to snow on land, land ice, and sea ice to produce consistent shortwave radiative properties for snow-covered surfaces across the Earth system. With the evolution and further understanding of snow and aerosol physics and chemistry, the adoption of this hybrid model will obviate the effort to modify and maintain separate optical variable input files used for different model components.

SNICAR-AD is now implemented in both the sea ice (MPAS-Seaice) and land (ELM) components of E3SM. More simulations and analyses are underway to examine its impact on E3SM model performance and simulated climate. The results are however beyond the scope of this work and will be thoroughly discussed in a future paper.

The data and models are available upon request to Cheng Dang (cdang5@uci.edu). SNICAR and dEdd–AD radiative transfer core can be found at https://github.com/E3SM-Project/E3SM (last access: 22 July 2019).

CD and CZ designed the study. CD coded the offline dEdd-AD and 2SD schemes, performed two-stream and 16-stream model simulations, and wrote the paper with input from CZ and MF. CZ performed the SWNB2 simulations. MF provided the base SNICAR code and snow optical inputs.

The authors declare that they have no conflict of interest.

This research has been supported by the U.S. Department of Energy (grant no. DE-SC0012998).

This paper was edited by Dirk Notz and reviewed by David Bailey and one anonymous referee.