Intercomparison and improvement of two-stream shortwave radiative transfer schemes in Earth system models for a unified treatment of cryospheric surfaces

Snow is an important climate regulator because it greatly increases the surface albedo of middle and high latitudes of the Earth. Earth system models (ESMs) often adopt two-stream approximations with different radiative transfer techniques, the same snow therefore has different solar radiative properties depending whether it is on land or on sea ice. Here we intercompare three two-stream algorithms widely used in snow models, improve their predictions at large zenith angles, and introduce a hybrid model suitable for all cryospheric surfaces in ESMs. The algorithms are those employed by the SNow ICe and Aerosol Radiative (SNICAR) module used in land models, dEdd–AD used in Icepack, the column physics used in the Los Alamos sea ice model CICE and MPAS-Seaice, and a two-stream discrete-ordinate (2SD) model. Compared with a 16-stream benchmark model, the errors in snow visible albedo for a direct-incident beam from all three two-stream models are small ($<\pm \mathrm{0.005}$) and increase as snow shallows, especially for aged snow. The errors in direct near-infrared (near-IR) albedo are small ($<\pm \mathrm{0.005}$) for solar zenith angles θ<75^∘, and increase as θ increases. For diffuse incidence under cloudy skies, dEdd–AD produces the most accurate snow albedo for both visible and near-IR ($<\pm \mathrm{0.0002}$) with the lowest underestimate (−0.01) for melting thin snow. SNICAR performs similarly to dEdd–AD for visible albedos, with a slightly larger underestimate (−0.02), while it overestimates the near-IR albedo by an order of magnitude more (up to 0.04). 2SD overestimates both visible and near-IR albedo by up to 0.03. We develop a new parameterization that adjusts the underestimated direct near-IR albedo and overestimated direct near-IR heating persistent across all two-stream models for θ>75^∘. These results are incorporated in a hybrid model SNICAR-AD, which can now serve as a unified solar radiative transfer model for snow in ESM land, land ice, and sea ice components.

1 Introduction

Snow cover on land, land ice, and sea ice, modulates the surface energy balance of middle and high latitudes of the Earth, principally because even a thin layer of snow can greatly increase the surface albedo. Integrated over the solar spectrum, the broadband albedo of opaque snow ranges from 0.7 to 0.9 (e.g., Wiscombe and Warren, 1980; Dang et al., 2015). In contrast, the albedo of other natural surfaces is smaller: 0.2, 0.25, and 0.5–0.7 for damp soil, grassland, and bare multi-year sea ice, respectively (Perovich, 1996; Liang et al., 2002; Brandt et al., 2005; Bøggild et al., 2010). The accumulation, evolution, and depletion of snow cover thus modify the seasonal cycle of surface albedo globally. In particular, snow over sea ice absorbs more solar energy and begins to melt in the spring, which forms melt ponds that bring the sea ice albedo to as low as 0.15 to further accelerate ice melt (Light et al., 2008, 2015). An accurate simulation of the shortwave radiative properties of snowpack is therefore crucial for spectrally partitioning solar energy and representing snow–albedo feedbacks across the Earth system. Unfortunately, computational demands and coupling architectures often constrain representation of snowpack radiative processes in Earth system models (ESMs; please refer to Table 1 for all abbreviations used in this work) to relatively crude approximations such as two-stream methods (Wiscombe and Warren, 1980; Toon et al., 1989). In this work, we intercompare two-stream methods widely used in snow models and then introduce a new parameterization that significantly reduces their snowpack reflectance and heating biases at large zenith angles, to produce more realistic behavior in polar regions.

Table 1 Abbreviations used in this paper and their references. Last access date for all cited URLs in this table is 22 July 2019.

Snow albedo is determined by many factors including the snow grain radius, the solar zenith angle, cloud transmittance, light-absorbing particles, and the albedo of underlying ground if snow is optically thin (Wiscombe and Warren, 1980; Warren and Wiscombe, 1980); it also varies strongly with wavelength since the ice absorption coefficient varies by 7 orders of magnitudes across the solar spectrum (Warren and Brandt, 2008). At visible wavelengths (0.2–0.7 µm), ice is almost nonabsorptive such that the absorption of visible energy by snowpack is mostly due to the light-absorbing particles (e.g., black carbon, organic carbon, mineral dust) that were incorporated during ice nucleation in clouds, scavenged during precipitation, or slowly sedimented from the atmosphere by gravity (Warren and Wiscombe, 1980, 1985; Doherty et al., 2010, 2014, 2016; Wang et al., 2013; Dang and Hegg, 2014). As snow becomes shallower, visible photons are more likely to penetrate through snowpack and get absorbed by darker underlying ground. At near-infrared (near-IR) wavelengths (0.7–5 µm), ice is much more absorptive, so that the snow near-IR albedo is lower than the visible albedo. Larger ice crystals form a lower albedo surface than smaller ice crystals; hence aged snowpacks absorb more solar energy. Photons incident at smaller solar zenith angles are more likely to penetrate deeper vertically and be scattered in the snowpack until being absorbed by the ice, the underlying ground, or absorbing impurities, which also leads to a smaller snow albedo. To compute the reflected solar flux, spectrally resolved albedo must be weighted by the incident solar flux, which is mostly determined by solar zenith angle, cloud cover and transmittance, and column water vapor. Modeling the solar properties of snowpacks must consider the spectral signatures of these atmospheric properties.

Several parameterizations have been developed to compute the snow solar properties without solving the radiative transfer equations and some are incorporated into ESMs or regional models. Marshall and Warren (1987) and Marshall (1989) parameterized snow albedo in both visible and near-IR bands as functions of snow grain size, solar zenith angle, cloud transmittance, snow depth, underlying surface albedo, and black carbon content. Marshall and Oglesby (1994) used this in an ESM. Gardner and Sharp (2010) computed the all-wave snow albedo with similar inputs. This was incorporated into the regional climate model RACMO (https://www.projects.science.uu.nl/iceclimate/models/racmo.php, last access: 22 July 2019) to simulate snow albedo in glaciered regions like Antarctica and Greenland (Kuipers Munneke et al., 2011). Dang et al. (2015) parameterized snow albedo as a function of snow grain radius, black carbon content, and dust content for visible and near-IR bands and 14 narrower bands used in the Rapid Radiative Transfer Model (RRTM; Mlawer and Clough, 1997). Their algorithm can also be expanded to different solar zenith angles using the zenith angle parameterization developed by Marshall and Warren (1987). Aoki et al. (2011) developed a more complex model based on the offline snow albedo and a transmittance look-up table. This can be applied to multilayer snowpack to compute the snow albedo and the solar heating profiles as functions of snow grain size, black carbon and dust content, snow temperature, and snowmelt water equivalent. These parameterizations are often in the form of simplified polynomial equations, which are especially suitable to long-term ESM simulations that require less time-consuming snow representations.

More complex models that explicitly solve the multiple-scattering radiative transfer equations have also been developed to compute snow solar properties. Flanner and Zender (2005) developed the SNow Ice and Aerosol Radiation model (SNICAR) that utilizes two-stream approximations (Wiscombe and Warren, 1980; Toon et al., 1989) to predict heating and reflectance for a multilayer snowpack. They implemented SNICAR in the Community Land Model (CLM) to predict snow albedo and vertically resolved solar absorption for snow-covered surfaces. Before SNICAR, CLM prescribed snow albedo and confined all solar absorption to the top snow layer (Flanner and Zender, 2005). Over the past decades, updates and new features have been added to SNICAR to consider more processes such as black carbon–ice mixing states (Flanner et al., 2012) and snow grain shape (He et al., 2018b). Concurrent with the development of SNICAR, Briegleb and Light (2007) improved the treatment of sea ice solar radiative calculations in the Community Climate System Model (CCSM). They implemented a different two-stream scheme with delta-Eddington approximation and the adding–doubling technique (hereafter, dEdd–AD) that allows CCSM to compute bare, ponded, and snow-covered sea ice albedo and solar absorption profiles of multilayer sea ice. Before these improvements, the sea ice albedo was computed based on surface temperature, snow thickness, and sea ice thickness using averaged sea ice and snow albedo. dEdd–AD has been adopted by the sea ice physics library Icepack (https://github.com/CICE-Consortium/Icepack/wiki, last access: 22 July 2019), which is used by the Los Alamos sea ice model CICE (Hunke et al., 2010) and Model for Prediction Across Scales Sea Ice (MPAS-Seaice; Turner et al., 2019). CICE itself is used in numerous global and regional models.

SNICAR and dEdd–AD solve the multiple-scattering radiative transfer equations and provide much improved solar radiative representations for the cryosphere, though their separate development and implementation created an artificial divide for snow simulation. In ESMs that utilize both SNICAR and dEdd–AD, such as the Community Earth System Model (CESM, http://www.cesm.ucar.edu/, last access: 22 July 2019) and the Energy Exascale Earth System Model (E3SM, previously known as ACME, https://e3sm.org/, last access: 22 July 2019), the solar radiative properties of snow on land and snow on sea ice are computed separately via SNICAR and dEdd–AD. As a result, the same snow in nature has different solar radiative properties such as reflectance depending on which model represents it. These differences are model artifacts that should be eliminated so that snow has consistent properties across the Earth system.

In this paper, we evaluate the accuracy and biases of three two-stream models listed in Table 2, including the algorithms used in SNICAR and dEdd–AD, for representing reflectance and heating. In Sects. 2–4, we describe the radiative transfer algorithms and calculations performed in this work. The results and model intercomparisons are discussed in Sect. 5. In Sect. 6, we introduce a parameterization to reduce the simulated albedo and heating bias for solar zenith angles larger than 75^∘. In Sect. 7, we summarize the major differences of algorithm implementations between SNICAR and dEdd–AD in ESMs. We use these results to develop and justify a unified surface shortwave radiative transfer method for all Earth system model components in the cryosphere, presented in Sect. 8.

Table 2 Two-stream radiative transfer algorithms evaluated in this work, including algorithms that are currently implemented in Earth system models CESM and E3SM.

2 Radiative transfer model

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.

2.1 SNICAR in land models CLM and ELM

SNICAR is implemented as the default snow shortwave radiative transfer scheme in CLM and the E3SM land model (ELM). It adopts the two-stream algorithms and the rapid solver developed by Toon et al. (1989) to compute the solar properties of multilayer snowpacks. These two-stream algorithms are derived from the general equation of radiative transfer in a plane-parallel media:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& \mathit{\mu }{\displaystyle \frac{\partial I}{\partial \mathit{\tau }}}\left(\mathit{\tau },\mathit{\mu },\mathrm{\Phi }\right)=I\left(\mathit{\tau },\mathit{\mu },\mathrm{\Phi }\right)-{\displaystyle \frac{\mathit{\varpi }}{\mathrm{4}\mathit{\pi }}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\int }_{-\mathrm{1}}^{\mathrm{1}}P\left(\mathit{\mu },\mathit{\mu }{}^{\prime },\mathit{\varphi },\mathit{\varphi }{}^{\prime }\right)\\ & I\left(\mathit{\tau },\mathit{\mu }{}^{\prime },\mathrm{\Phi }{}^{\prime }\right)\mathrm{d}\mathit{\mu }{}^{\prime }\mathrm{d}\mathit{\varphi }{}^{\prime }-S(\mathit{\tau },\mathit{\mu },\mathrm{\Phi })\end{array},\end{array}$$

where Φ is azimuth angle, μ is the cosine of the zenith angle, and ϖ is single-scattering albedo. On the right-hand side, the three terms are intensity at optical depth τ, internal source term due to multiple scattering, and external source term S. For a purely external source at solar wavelengths S is

$$\begin{array}{}\text{(2)}& S={\displaystyle \frac{\mathit{\varpi }}{\mathrm{4}}}{F}_{\mathrm{s}}P\left(\mathit{\mu },-{\mathit{\mu }}_{\mathrm{0}},\mathit{\varphi },{\mathit{\varphi }}_{\mathrm{0}}\right)\exp \left({\displaystyle \frac{-\mathit{\tau }}{{\mathit{\mu }}_{\mathrm{0}}}}\right),\end{array}$$

where πF_s is incident solar flux, and μ₀ is the incident direction of the solar beam. Integrating Eq. (1) over azimuth and zenith angles yields the general solution of two-stream approximations (Meador and Weaver, 1980). The upward and downward fluxes at optical depth τ of layer n can be represented as

$$\begin{array}{}\text{(3a)}& {\displaystyle }& {\displaystyle }{F}_n^{+}=k_{\mathrm{1}n}\exp \left({\mathrm{\Lambda }}_n\mathit{\tau }\right)+{\mathrm{\Gamma }}_n{k}_{\mathrm{2}n}\exp \left(-{\mathrm{\Lambda }}_n\mathit{\tau }\right)+C_n^{+}\left(\mathit{\tau }\right),\text{(3b)}& {\displaystyle }& {\displaystyle }{F}_n^{-}={\mathrm{\Gamma }}_n{k}_{\mathrm{1}n}\exp \left({\mathrm{\Lambda }}_n\mathit{\tau }\right)+k_{\mathrm{2}n}\exp \left(-{\mathrm{\Lambda }}_n\mathit{\tau }\right)+C_n^{-}\left(\mathit{\tau }\right),\end{array}$$

where Λ_n, Γ_n, and C_n are known coefficients determined by the two-stream method, incident solar flux, and solar zenith angle; whereas k_1n and k_2n are unknown coefficients determined by the boundary conditions. For an N-layer snowpack, the solutions for upward and downward fluxes are coupled at layer interfaces to generate 2N equations with 2N unknown coefficients k_1n and k_2n. Combining these equations linearly generates a new set of equations with terms in tri-diagonal form that enables the application of a fast tri-diagonal matrix solver. With the solved coefficients, the upward and downward fluxes are computed at different optical depths (Eqs. 3a and 3b) and eventually the reflectance, transmittance, and absorption profiles of solar flux for any multilayer snowpack.

SNICAR itself implements all three two-stream algorithms in Toon et al. (1989): Eddington, quadrature, and hemispheric mean. In practical simulations, it utilizes the Eddington and hemispheric-mean approximations to compute the visible and near-IR snow properties, respectively (Flanner et al., 2007). In addition to its algorithms, SNICAR implements the delta transform of the fundamental input variable asymmetry factor (g), single-scattering albedo (ϖ), and optical depth (τ) to account for the strong forward scattering in snow (Eqs. 2a–2c, Wiscombe and Warren, 1980).

2.2 dEdd–AD in sea ice models Icepack, CICE, and MPAS-Seaice

Icepack, CICE, and MPAS-Seaice use the same shortwave radiative scheme dEdd–AD developed and documented by Briegleb and Light (2007). Sea ice is divided into multiple layers to first compute the single-layer reflectance and transmittance using two-stream delta-Eddington solutions to account for the multiple scattering of light within each layer (Equation set 50, Briegleb and Light, 2007), where the name “delta” implies dEdd–AD implements the delta transform to account for the strong forward scattering of snow and sea ice (Eqs. 2a–2c, Wiscombe and Warren, 1980). The single-layer direct albedo and transmittance are computed by equations

$$\begin{array}{}\text{(4a)}& {\displaystyle }\begin{array}{rl}& R\left({\mathit{\mu }}_{\mathrm{0},n}\right)=A_n\exp \left({\displaystyle \frac{-\mathit{\tau }}{{\mathit{\mu }}_{\mathrm{0},n}}}\right)\\ & +B_n\left(\exp \left({\mathit{\epsilon }}_n\mathit{\tau }\right)-\exp \left(-{\mathit{\epsilon }}_n\mathit{\tau }\right)\right)-K_n,\end{array}\text{(4b)}& {\displaystyle }\begin{array}{rl}& T\left({\mathit{\mu }}_{\mathrm{0},n}\right)=E_n\\ & +H_n\left(\exp \left({\mathit{\epsilon }}_n\mathit{\tau }\right)-\exp \left(-{\mathit{\epsilon }}_n\mathit{\tau }\right)\right)\exp \left({\displaystyle \frac{-\mathit{\tau }}{{\mathit{\mu }}_{\mathrm{0},n}}}\right),\end{array}\end{array}$$

where coefficients A_n, B_n, K_n, E_n, H_n, and ε_n are determined by the single-scattering albedo (ϖ), asymmetry factor (g), optical depth (τ), and angle of the incident beam at layer n (μ_0, n). Following the delta-Eddington assumption, simple formulas are available for the single-layer reflectance and transmittance under both clear sky (direct flux, Eqs. 4a and 4b) and overcast sky (diffuse flux) conditions. However, the formula derived by applying diffuse-flux upper boundary conditions sometimes yields negative albedos (Wiscombe, 1977). To avoid the unphysical values, diffuse reflectance $\stackrel{\mathrm{‾}}{R}$ and transmittance $\stackrel{\mathrm{‾}}{T}$ of a single layer are computed by integrating the direct reflectance R(μ) and transmittance T(μ) over the incident hemisphere assuming isotropic incidence:

$$\begin{array}{}\text{(5a)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{R}=\mathrm{2}{\int }_{\mathrm{0}}^{\mathrm{1}}\mathit{\mu }R\left(\mathit{\mu }\right)\mathrm{d}\mathit{\mu },\text{(5b)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{T}=\mathrm{2}{\int }_{\mathrm{0}}^{\mathrm{1}}\mathit{\mu }T\left(\mathit{\mu }\right)\mathrm{d}\mathit{\mu }.\end{array}$$

This is the same as the method proposed by Wiscombe and Warren (1980, their Eq. 5). In practice, eight Gaussian angles are implemented to perform the integration for every layer.

The computed single-layer reflectance and transmittance of direct and diffuse components are then combined to account for the interlayer scattering of light to compute the reflectance and transmission at every interface (Equation set 51, Briegleb and Light, 2007), and eventually the upward and downward fluxes (Equation set 52, Briegleb and Light, 2007). These upward and downward fluxes at each optical depth are then used to compute the column reflectance and transmittance, and the absorption profiles for any multilayered media, such as snowpacks on land and sea ice.

In nature, a large fraction of sea ice is covered by snow during winter. As snow melts away in late spring and summer, it exposes bare ice, and melt ponds form on the ice surface. Such variation in sea ice surface types requires the shortwave radiative transfer model to be flexible and capable of capturing the light refraction and reflection. Refractive boundaries exist where air (refractive index m_re=1.0), snow (assuming snow as medium of air containing a collection of ice particles, m_re=1.0), pond (assuming pure water, m_re=1.33), and ice (assuming pure ice, m_re=1.31) are present in the same sea ice column. The general solution of delta-Eddington and the two-stream algorithms used in SNICAR are not applicable to such nonuniformly refractive layered media. To include the effects of refraction, Briegleb and Light (2007) modified the adding formula at the refractive boundaries (i.e., interfaces between air and ice, snow and ice, and air and pond). The reflectance and transmittance of the adjacent layers above and below the refractive boundary are combined with modifications to include the Fresnel reflection and refraction of direct and diffuse fluxes (Sect. 4.1, Briegleb and Light, 2007). dEdd–AD can thus be applied to any layered media with either uniform (e.g., snow on land) or nonuniform (e.g., snow on sea ice) refractive indexes.

In this paper, we apply dEdd–AD to snowpacks that can be treated as uniform refractive media such as the land snow columns assumed in SNICAR for model evaluation. An ideal radiative treatment for snow should, however, keep the potential to include refraction for further applications to snow on sea ice or ice sheets. Therefore, in addition to these two widely used algorithms in Icepack and SNICAR, we evaluate a third algorithm (Sect. 2.3) that can be applied to layered media with either uniform or nonuniform refractive indexes.

2.3 Two-stream discrete-ordinate algorithm (2SD)

A refractive boundary also exists between the atmosphere and the ocean, and models have been developed to solve the radiative transfer problems in the atmosphere–ocean system using the discrete-ordinate technique (e.g., Jin and Stamnes, 1994; Lee and Liou, 2007). Similar to the two-stream algorithms of Toon et al. (1989) used in SNICAR, Jin and Stamnes (1994) also developed their algorithm from the general equation

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& \mathit{\mu }{\displaystyle \frac{\partial I}{\partial \mathit{\tau }}}\left(\mathit{\tau },\mathit{\mu }\right)=I\left(\mathit{\tau },\mathit{\mu }\right)\\ & -{\displaystyle \frac{\mathit{\varpi }}{\mathrm{4}\mathit{\pi }}}{\int }_{-\mathrm{1}}^{\mathrm{1}}P\left(\mathit{\tau },\mathit{\mu },\mathit{\mu }{}^{\prime }\right)I\left(\mathit{\tau },\mathit{\mu }{}^{\prime }\right)\mathrm{d}\mathit{\mu }{}^{\prime }-S\left(\mathit{\tau }\mathit{\mu }\right)\end{array}.\end{array}$$

Equation (6) is the azimuthally integrated version of Eq. (1). However, for vertically inhomogeneous media like the atmosphere–ocean or sea ice, the external source term S(τ, μ) is different. Specifically, for the medium of total optical depth τ^a above the refractive interface, one must consider the contribution from the upward beam reflected at the refractive boundary (second term on the right-hand side):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& S^{\mathrm{a}}\left(\mathit{\tau },\mathit{\mu }\right)={\displaystyle \frac{\mathit{\varpi }}{\mathrm{4}\mathit{\pi }}}{F}_{\mathrm{s}}P\left(\mathit{\tau },-{\mathit{\mu }}_{\mathrm{0}},\mathit{\mu }\right)\exp \left({\displaystyle \frac{-\mathit{\tau }}{{\mathit{\mu }}_{\mathrm{0}}}}\right)\\ & +{\displaystyle \frac{\mathit{\varpi }}{\mathrm{4}\mathit{\pi }}}{F}_{\mathrm{s}}R\left(-{\mathit{\mu }}_{\mathrm{0}},m\right)P\left(\mathit{\tau },+{\mathit{\mu }}_{\mathrm{0}},\mathit{\mu }\right)\exp \left({\displaystyle \frac{-\left(\mathrm{2}{\mathit{\tau }}^{\mathrm{a}}-\mathit{\tau }\right)}{{\mathit{\mu }}_{\mathrm{0}}}}\right)\end{array},\end{array}$$

where $R\left(-{\mathit{\mu }}_{\mathrm{0}},m\right)$ is the Fresnel reflectance of radiation and m is the ratio of the refractive indices of the lower to the upper medium. For the medium below the refractive interface, one must account for the Fresnel transmittance $T\left(-{\mathit{\mu }}_{\mathrm{0}},m\right)$ and modify the angle of beam travel in media b:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& S^{\mathrm{b}}\left(\mathit{\tau },\mathit{\mu }\right)={\displaystyle \frac{\mathit{\varpi }}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{{\mathit{\mu }}_{\mathrm{0}n}}}{F}_{\mathrm{s}}T\left(-{\mathit{\mu }}_{\mathrm{0}},m\right)P\left(\mathit{\tau },-{\mathit{\mu }}_{\mathrm{0}},\mathit{\mu }\right)\\ & \exp \left({\displaystyle \frac{-{\mathit{\tau }}^{\mathrm{a}}}{{\mathit{\mu }}_{\mathrm{0}}}}\right)\exp \left({\displaystyle \frac{-\left(\mathit{\tau }-{\mathit{\tau }}^{\mathrm{a}}\right)}{{\mathit{\mu }}_{\mathrm{0}n}}}\right)\end{array},\end{array}$$

where μ_0n is the cosine zenith angle of refracted beam incident at angle μ₀ above the refractive boundary, by Snell's law:

$$\begin{array}{}\text{(9)}& {\mathit{\mu }}_{\mathrm{0}n}=\sqrt{\mathrm{1}-\left(\mathrm{1}-{\mathit{\mu }}_{\mathrm{0}}^{\mathrm{2}}\right)/m^{\mathrm{2}}}.\end{array}$$

For uniformly refractive media like snow on land, one can just set the refractive index m_re equal to 1 for every layer. In this case, the Fresnel reflectance $R\left(-{\mathit{\mu }}_{\mathrm{0}},m\right)$ is 0 in Eq. (7), the Fresnel transmittance $T\left(-{\mathit{\mu }}_{\mathrm{0}},m\right)$ is 1 in Eq. (8), and μ_0n equals μ₀: the two source terms S^a(τ, μ) and S^b(τ, μ) become the same and equal the source term of homogenous media given in Eq. (2).

For two-stream approximations of this method, analytical solutions of upward and downward fluxes are coupled at each layer interface to generate 2N equations with 2N unknown coefficients for any N-layer stratified column. The solutions of two-stream algorithms and boundary conditions for homogenous media are well documented (Sect. 8.4 and 8.10 of Thomas and Stamnes, 1999). Despite the extra source terms, these 2N equations can also be organized into a tri-diagonal matrix similar to the method of Toon et al. (1989) used in SNICAR. Flexibility and speed therefore make this two-stream discrete-ordinate algorithm (hereafter, 2SD) a potentially good candidate for long-term Earth system modeling. In this work, we only apply 2SD to the snowpack and note that it can be applied to any uniformly or nonuniformly refractive media like snow on land or sea ice, with the delta transform implemented for fundamental optical variables (Eqs. 2a–2c, Wiscombe and Warren, 1980).

2.4 16-stream DISORT

In addition to the mathematical technique, the accuracy and speed of radiative transfer algorithms depend on the number of angles used for flux estimation in the upward and downward hemispheres. SNICAR, dEdd–AD, and 2SD use one angle to represent upward flux and one angle to represent downward flux; hence they are named the two-stream algorithm. Lee and Liou (2007) use two upward and two downward streams. Jin and Stamnes (1994) documented the solutions for any even number of streams. The computational efficiency of these models is lower than that of two-stream models while their accuracy is better. To quantify the accuracy of the three two-stream algorithms for snow shortwave simulations, we use the 16-stream DIScrete-Ordinate Radiative Transfer model (DISORT) as the benchmark model (http://lllab.phy.stevens.edu/disort/, last access: ) (Stamnes et al., 1988).

3 Input for radiative transfer models

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:

$$\begin{array}{}\text{(10)}& \mathit{\beta }=\mathrm{3}/\left(r{\mathit{\rho }}_{\mathrm{ice}}\right),\end{array}$$

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.

4 Solar spectra used for the spectral integrations

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 α(λ):

$$\begin{array}{}\text{(11)}& \mathit{\alpha }={\displaystyle \frac{{\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}\mathit{\alpha }\left(\mathit{\lambda }\right)F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{{\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}}.\end{array}$$

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 1 Spectral 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.

5 Model evaluation

5.1 Spectral albedo and reflected solar flux

The spectral reflectance of pure deep snow computed using two-stream models and 16-stream DISORT is shown in Fig. 2. The snow grain radius is 100 µm – a typical grain size for fresh new snow. For clear sky with a direct beam source (left column), all three two-stream models show good accuracy at visible wavelengths (0.3–0.7 µm), and within this band, the snow albedo is large and close to 1. As wavelength increases, the albedo diminishes in the near-IR band. Two-stream models overestimate snow albedo at these wavelengths, with maximum biases of 0.013 (SNICAR and dEdd–AD) and 0.023 (2SD) within wavelength 1–1.7 µm. For cloudy-sky cases with diffuse upper boundary conditions, dEdd–AD reproduces the snow albedo at all wavelengths with the smallest absolute error (<0.005), and SNICAR and 2SD both overestimate the snow albedo with maximum biases >0.04 between 1.1 and 1.4 µm.

Figure 2 The spectral albedo of pure snow computed using 16-stream DISORT, SNICAR, dEdd–AD, and 2SD models, for clear-sky (direct beam at solar zenith angle 60^∘) and cloudy-sky conditions in the left and right panels, respectively. Panels (a, b) show spectral albedo. Panels (c, d) show the difference ($\mathit{\delta }\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{2}}-{\mathit{\alpha }}_{\mathrm{16}}$) in spectral albedos computed using the two-stream model (α₂) and 16-stream DISORT (α₁₆). Panels (e, f) show the difference of reflected spectral flux given δα. The snowpack is set to semi-infinite deep with a grain radius of 100 µm.

In both sky conditions, the errors of snow albedo are larger at near-IR wavelengths ranging from 1.0 to 1.7 µm, while the solar incident flux peaks at 0.5 µm then decrease as wavelength increases. The largest error in reflected flux is within the 0.7–1.5 µm band for SNICAR and 2SD, as shown in the third row of Fig. 2. dEdd–AD overestimates the direct snow albedo mostly at wavelengths larger than 1.5 µm where the error in reflected flux is almost negligible.

5.2 Broadband albedo and reflected solar flux

Integrated over the visible and near-IR wavelengths, the error in band albedos computed using two-stream models for different cases is shown in Figs. 3–6.

Figure 3 The difference in direct snow albedo ($\mathit{\delta }\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{2}}-{\mathit{\alpha }}_{\mathrm{16}}$) computed using two-stream models (α₂) and using the 16-stream DISORT model (α₁₆), for various snow depths and solar zenith angles, with a snow grain radius of 100 µm. From the top to the bottom, rows are results of two-stream models SNICAR, dEdd–AD, and 2SD. From the left to the right columns are albedo differences of all-wave, visible, and near-IR bands.

Figure 3 shows the error in direct band albedo for fixed snow grain radius of 100 µm with different snow depth and solar zenith angles. As introduced in Sect. 2, SNICAR and dEdd–AD both use the delta-Eddington method to compute the visible albedo. They overestimate the visible albedo for solar zenith angles smaller than 50^∘ by up to 0.005, and underestimate it for solar zenith angles larger than 50^∘ by up to −0.01. 2SD produces similar results for the visible band but at a larger solar zenith angle threshold of 75^∘. In the near-IR band, SNICAR and 2SD overestimate the snow albedo for solar zenith angles smaller than 70^∘, beyond this, the error in albedo increases by up to −0.1 as solar zenith angle increases. dEdd–AD produces a similar error pattern with a smaller solar zenith angle threshold at 60^∘. As snow ages, its average grain size increases. For typical old melting snow of grain radius 1000 µm (Fig. 4), two-stream models produce similar errors of direct albedo in all bands. Integrating over the entire solar band, the three two-stream models evaluated show similar error patterns for direct albedo.

Figure 4 The difference in direct snow albedo ($\mathit{\delta }\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{2}}-{\mathit{\alpha }}_{\mathrm{16}}$) computed using two-stream models (α₂) and using the 16-stream DISORT model (α₁₆), for various snow depths and solar zenith angles, with a snow grain radius of 1000 µm.

For a fixed solar zenith angle of 60^∘, the error of direct albedo for different snow depth and snow grain radii is shown in Fig. 5. SNICAR and dEdd–AD underestimate the visible albedo in most scenarios, while 2SD overestimates the visible albedo for a larger range of grain radius and snow depth. All three two-stream models tend to overestimate the near-IR albedo except for shallow snow with large grain radius; the error of 2SD is 1 order of magnitude larger than that of SNICAR and dEdd–AD.

Figure 5 The difference in direct snow albedo ($\mathit{\delta }\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{2}}-{\mathit{\alpha }}_{\mathrm{16}}$) computed using two-stream models (α₂) and using the 16-stream DISORT model (α₁₆), for various snow depths and snow grain radii, with a solar zenith angle of 60^∘.

Figure 6 is similar to Fig. 5, but shows the diffuse snow albedo. In the visible band, SNICAR and dEdd–AD generate similar errors in that they both underestimate the albedo as snow grain size increases and snow depth decreases. 2SD overestimates the albedo with a maximum error of around 0.015. In the near-IR, two-stream models tend to overestimate snow albedo, while the magnitude of biases produced by SNICAR and 2SD is 1 order larger than that of dEdd–AD with the maximum error of 0.035 generated by SNICAR. As a result, the all-wave diffuse albedos computed using dEdd–AD are more accurate than those computed using SNICAR and 2SD.

Figure 6 The difference in diffuse snow albedo ($\mathit{\delta }\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{2}}-{\mathit{\alpha }}_{\mathrm{16}}$) computed using two-stream models (α₂) and using the 16-stream DISORT model (α₁₆), for various snow depths and snow grain radii, with a solar zenith angle of 60^∘ at the top of the atmosphere.

Figures 7, 8, and 9 show the errors in reflected shortwave flux caused by snow albedo errors seen in Figs. 3, 4, and 6. In general, two-stream models produce larger errors in reflected direct near-IR flux (Figs. 7 and 8), especially with the 2SD model: the maximum overestimate of reflected near-IR flux is 6–8 W m⁻² for deep melting snow with a solar zenith angle <30^∘. Errors in reflected direct visible flux are smaller (mostly within ±1 W m⁻²) for all models in most scenarios, and become larger (mostly within ±3 W m⁻²) as snow grain size increases to 1000 µm if computed using 2SD. As shown in Fig. 9, for diffuse flux with a solar zenith angle of 60^∘ at the top of the atmosphere (TOA), SNICAR and dEdd–AD generate small errors in reflected visible flux (mostly within ±1 W m⁻²), while 2SD always overestimates reflected visible flux by up to 5 W m⁻². In the near-IR, SNICAR and 2SD overestimate reflected flux by as much as 10–12 W m⁻²; the error in reflected near-IR flux produced by dEdd–AD is much smaller, mostly within ±1 W m⁻².

Figure 7 Error in reflected direct solar flux given albedo errors shown in Fig. 3.

Figure 8 Error in reflected direct solar flux given albedo errors shown in Fig. 4.

Figure 9 Error in reflected diffuse solar flux given albedo errors shown in Fig. 6.

In general, dEdd–AD produces the most accurate albedo and thus reflected flux for both direct and diffuse components. SNICAR is similar to dEdd–AD for its accuracy of direct albedo and flux, yet generates large error for the diffuse component. 2SD tends to overestimate snow albedo and reflected flux in both direct and diffuse components and shows the largest errors among three two-stream models. Although the differences between algorithms are small, they can have a notable impact on snowpack melt. For example, compared to dEdd–AD, SNICAR and 2SD overestimate the diffuse albedo by ∼0.015 for melting snow (Fig. 6). In Greenland, the daily averaged downward diffuse solar flux from May to September is 200 W m⁻², and the averaged cloud cover fraction is 80 % (Fig. 6, Dang et al., 2017). In this case, SNICAR and 2SD overestimate the reflected solar flux by 2.4 W m⁻² d⁻¹ – the amount of energy is otherwise enough to melt 10 cm of snow water equivalent from May to September. dEdd–AD also remediates compensating spectral biases (where visible and near-IR biases are of opposite signs) present in the other schemes. Those spectral biases do not affect the broadband fluxes like the diffuse biases, but they nevertheless degrade proper feedbacks between snow–ice reflectance and heating.

5.3 Band absorption of solar flux

Figure 10 shows absorption profiles of shortwave flux computed using the 16-stream DISORT model, with errors in absorbed fractional solar flux computed using two-stream models. The snowpack is 10 cm deep and is divided into five layers, each 2 cm thick. The snow grain radii are set to 100 µm and 1000 µm. The figure shows fractional absorption for snow layers 1–4 and the underlying ground with an albedo of 0.25.

As shown in the first column of Fig. 10, for new snow with a radius of 100 µm, most solar absorption occurs in the top 2 cm snow layer, where roughly 10 % and 15 % of diffuse and direct near-IR flux is absorbed and dominates the solar absorption within the snowpack. In the second layer (2–4 cm), the absorption of solar flux is less than 1 % and gradually decreases within the interior layers. The underlying ground absorbs roughly 2 % of solar flux, mostly visible flux that penetrates the snowpack more efficiently. As snow ages and snow grain grows, photons penetrate deeper into the snowpack. For typical old melting snow with a radius of 1000 µm, most solar absorption still occurs in the top 2 cm snow layer, where roughly 20 % and 14 % of diffuse and direct near-IR flux is absorbed. The second snow layer (2–4 cm) absorbs more near-IR solar flux by roughly 2 %. More photons can penetrate through the snowpack, and result in a high fractional absorption by the underlying ground, especially for the visible band. As snow depth increases, the ground absorption will decrease for both snow radii.

Figure 10 Comparison of light-absorption profiles derived from two-stream models and 16-stream DISORT. The left-most column shows fractional band absorptions computed using 16-stream DISORT. The right three panels show the errors of all-wave, visible, and near-IR fractional absorptions calculated using two-stream models. The top and bottom panels are for clear-sky and cloudy-sky conditions (solar zenith angle of 60^∘), respectively. The snowpack is 10 cm deep and is divided evenly into five 2 cm thick layers, for new snow (r=100 µm) and old snow (r=1000 µm). Layers 1–4 represent the top four snow layers (top 8 cm), and layer 5 represents underlying ground with an albedo of 0.25.

Comparing to 16-stream DISORT, two-stream models underestimate the column solar absorptions for new snow, and they overestimate them for old snow, especially for the surface snow layer and the underground. Overall, dEdd–AD gives the most accurate absorption profiles among the three two-stream models, especially for new snow.

6 Correction for direct albedo for large solar zenith angles

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.

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:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& R_{{\mathrm{75}}_{+}}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{16}\text{-}\mathrm{DISORT}}}{{\mathit{\alpha }}_{\mathrm{dEdd}-\mathrm{AD}}}}=c_{\mathrm{1}}\left({\mathit{\mu }}_{\mathrm{0}}\right){\log }_{\mathrm{10}}\left(r\right)+c_{\mathrm{0}}\left({\mathit{\mu }}_{\mathrm{0}}\right),\\ & \text{for}{\mathit{\mu }}_{\mathrm{0}}<\mathrm{0.26},\text{i.e.,}{\mathit{\theta }}_{\mathrm{0}}>\mathrm{75}{}^{\circ },\end{array}\end{array}$$

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

$$\begin{array}{}\text{(13a)}& {\displaystyle }& {\displaystyle }{c}_{\mathrm{1}}\left({\mathit{\mu }}_{\mathrm{0}}\right)=\mathrm{1.304}{\mathit{\mu }}_{\mathrm{0}}^{\mathrm{2}}-\mathrm{0.631}{\mathit{\mu }}_{\mathrm{0}}+\mathrm{0.086},\text{(13b)}& {\displaystyle }& {\displaystyle }{c}_{\mathrm{0}}\left({\mathit{\mu }}_{\mathrm{0}}\right)=\mathrm{6.807}{\mathit{\mu }}_{\mathrm{0}}^{\mathrm{2}}-\mathrm{3.338}{\mathit{\mu }}_{\mathrm{0}}+\mathrm{1.467}.\end{array}$$

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}}_{+}}$:

$$\begin{array}{}\text{(14a)}& {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\mathrm{dEdd}\text{-}\mathrm{AD}}^{\mathrm{adjust}}=R_{{\mathrm{75}}_{+}}{\mathit{\alpha }}_{\mathrm{dEdd}\text{-}\mathrm{AD}},\text{(14b)}& {\displaystyle }& {\displaystyle }{\mathrm{\text{Fabs}}}_{\mathrm{dEdd}\text{-}\mathrm{AD}}^{\mathrm{adjust}}={\text{Fabs}}_{\mathrm{dEdd}\text{-}\mathrm{AD}}-(R_{{\mathrm{75}}_{+}}-\mathrm{1}){\mathit{\alpha }}_{\mathrm{dEdd}\text{-}\mathrm{AD}}{F}_{\mathrm{nir}},\end{array}$$

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 12 Error 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.

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.

7 Implementation of snow radiative transfer model in Earth system models

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:

$$\begin{array}{}\text{(15a)}& {\displaystyle }& {\displaystyle }\mathit{\varpi }\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)={\displaystyle \frac{{\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}\mathit{\varpi }\left(\mathit{\lambda }\right)F\left(\mathit{\lambda }\right)\mathit{\alpha }\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{{\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}F\left(\mathit{\lambda }\right)\mathit{\alpha }\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},\text{(15b)}& {\displaystyle }& {\displaystyle }g\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)={\displaystyle \frac{{\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}g\left(\mathit{\lambda }\right)F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{{\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}F\left(\mathit{\lambda }\right)\mathit{\alpha }\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},\text{(15c)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{ext}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)={\displaystyle \frac{{\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}{\mathit{\sigma }}_{\mathrm{ext}}\left(\mathit{\lambda }\right)F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{{\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}F\left(\mathit{\lambda }\right)\mathit{\alpha }\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}}.\end{array}$$

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:

$$\begin{array}{}\text{(16a)}& {\displaystyle }& {\displaystyle }{w}_{\mathrm{nir}\mathrm{1}}=\mathrm{0.67}+\mathrm{0.11}\cdot \left(\mathrm{1}-f_{\mathrm{nidr}}\right),\text{(16b)}& {\displaystyle }& {\displaystyle }{w}_{\mathrm{nir}\mathrm{2}}=\mathrm{1}-w_{\mathrm{nir}\mathrm{1}}.\end{array}$$

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.

$$\begin{array}{}\text{(17a)}& {\displaystyle }& {\displaystyle }\mathit{\varpi }\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)={\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}\mathit{\varpi }\left(\mathit{\lambda }\right)F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\text{(17b)}& {\displaystyle }& {\displaystyle }g\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)={\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}g\left(\mathit{\lambda }\right)F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\text{(17c)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{ext}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)={\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}{\mathit{\sigma }}_{\mathrm{ext}}\left(\mathit{\lambda }\right)F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\end{array}$$

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

$$\begin{array}{}\text{(18)}& \stackrel{\mathrm{‾}}{\mathit{\alpha }}=\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}\mathit{\alpha }\left(\mathit{\lambda }\right)F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }.\end{array}$$

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):

$$\begin{array}{}\text{(19)}& r_{\mathrm{eqv}}=\left(f_{\mathrm{nidr}}+\mathrm{0.8}\left(\mathrm{1}-f_{\mathrm{nidr}}\right)\right)r.\end{array}$$

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.

8 Towards a unified radiative transfer model for snow, sea ice, and land ice

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.

9 Conclusions

In this work, we aim to improve and unify the solar radiative transfer calculations for snow on land and snow on sea ice in ESMs by evaluating the following two-stream radiative transfer algorithms: the two-stream delta-Eddington adding–doubling algorithm dEdd–AD implemented in sea ice models Icepack, CICE, and MPAS-Seaice, the two-stream delta-Eddington and two-stream delta-Hemispheric-Mean algorithms implemented in snow model SNICAR, and a two-stream delta-discrete-ordinate algorithm. Among these three models, dEdd–AD produces the most accurate snow albedo and solar absorption (Sect. 5). All two-stream models underestimate near-IR snow albedo and overestimate near-IR absorption when solar zenith angles are larger than 75^∘, which can be adjusted by a parameterization we developed (Sect. 6). We compared the implementations of radiative transfer cores in SNICAR and dEdd–AD (Sect. 7) and recommended a consistent and hybrid shortwave radiative model SNICAR-AD for snow-covered surfaces across ESMs (Sect. 8). Improved treatment of surface cryospheric radiative properties in the thermal infrared has recently been shown to remediate significant climate simulation biases in polar regions (Huang et al., 2018). It is hoped that adoption of improved and consistent treatments of solar radiative properties for snow-covered surfaces as described in this study will further remediate simulation biases in snow-covered regions.