Improvement of a Snow Albedo Parameterization in the Snow-Atmosphere-Soil Transfer Model: Evaluation of Impacts of Aerosol on Seasonal Snow Cover

Snow is an important component in the cryosphere, which plays a significant role in the global climate system by exerting effects on heat absorption by the surface and providing the water available for evaporation into the atmosphere (e.g., Robock, 1980; Barnett et al., 1989; Groisman et al., 1994; Sun et al., 1999; Oaida et al., 2015). A number of snow-process schemes in land-surface and climate models have been developed to describe the mass and energy balance of the snowpack and the interaction between the surface and the atmosphere (e.g., Anderson, 1976; Jordan, 1991; Verseghy, 1991; Loth et al., 1993; Sun et al., 1999; Oleson et al., 2010; Vionnet et al., 2012). The Inter-comparison of Land Surface Parameterization Schemes Project compared different land-surface schemes used for climate and weather forecast models (Henderson-Sellers et al., 1993). It was found that most of them capture the basic features of snowmelt and runoff reasonably. However, land-surface models with different levels of complexity exhibit large discrepancies in their simulations of snow accumulation, snow melting, turbulent heat fluxes, and streamflow (Henderson-Sellers et al., 1993; Nijssen et al., 2003). These discrepancies may be attributable to the differences in various internal snow parameterizations and interactions among the model components.

One such parameterization, which critically modulates the energy balance between the surface and the atmosphere, is the snow-albedo parameterization. In addition, snow-albedo feedback is one of the most important climate feedbacks that could exert pronounced influences on the climate (e.g., Robock, 1983; Randall et al., 1994; Yang et al., 2001; Bony et al., 2006; Qu and Hall, 2007). (Thackeray and Fletcher, 2016) found that 1^°C of temperature rise roughly corresponds to about a 1% decrease in surface albedo when the feedback reaches its peak during springtime over the extratropics of the Northern Hemisphere. Apart from the climatic feedback, changes in the snow albedo could influence the snowmelt regime, which has major implications for water resources and ecological systems (e.g., Barnett et al., 2005; Steltzer et al., 2009; Painter et al., 2010; Bryant et al., 2013).

(Wiscombe and Warren, 1980) showed that snow albedo is related to the sky conditions and the physical properties of snow itself (i.e., snow grain size, liquid water content, snow depth, and impurity contents) (Warren and Wiscombe, 1980). In particular, changes in snow grain size influence the variation in the near-infrared albedo. Old and large snow grains will increase the path length of radiation passing through it and result in more radiation being absorbed by the snow, which leads to lower snow albedo. In addition, studies have reported that light-absorbing aerosols (LAA), particularly the light-absorbing dusts, black carbon (BC) and organic carbon deposited on the snowpack, are the dominant absorbers of radiation in the visible spectrum (e.g., Hansen and Nazarenko, 2004; McConnell et al., 2007; Ramanathan and Carmichael, 2008; Flanner et al., 2009; Skiles et al., 2012).

With the increasing emissions of aerosols due to human activity, increasing attention has been paid in recent decades on investigating the radiative forcing of LAA in snow (e.g., Hansen and Nazarenko, 2004; Randerson et al., 2006; Flanner et al., 2007; Oaida et al., 2015; Qian et al., 2015). (Bond et al., 2013) showed that the global effective BC forcing in snow and ice was +0.13 W m^-2. (Di Mauro et al., 2015) reported that the instantaneous radiative forcing of mineral dust was up to +153 W m^-2 at the most concentrated sampling area. Through such radiative forcing, the LAA in snow have a large impact on snow melting, snow duration, and runoff (e.g., Flanner and Zender, 2005; Painter et al., 2007, 2010; Flanner et al., 2009; Qian et al., 2009). (Painter et al., 2010) showed that peak runoff of the Colorado River at Lees Ferry, Arizona, occurred on average three weeks earlier than usual under heavier dust loading. Furthermore, the LAA in snow have major implications for the hydrologic cycle and atmospheric processes at different temporal and spatial scales in multiple ways (e.g., Xu et al., 2009; Qian et al., 2011, 2015; Bond et al., 2013; Oaida et al., 2015). (Xu et al., 2009) found that aerosols deposited on glaciers and snow cover contribute significantly to the observed rapid glacial retreat over the Tibetan Plateau (TP). LAA have been suggested to play an important role in TP climate change via affecting the heating of the snowpack, which can further influence the South and East Asian monsoons (Qian et al., 2011).

However, most of snow-albedo parameterizations employed in current land-surface and climate models neglect the effects of LAA in the snowpack. These schemes either just keep the snow albedo to a fixed value, or only make it an empirical function of snow depth, temperature, or snow age (e.g., Jordan, 1991; Loth et al., 1993; Sun et al., 1999; Roeckner et al., 2003), which is often the largest source of error in the simulation of the mass balance and snow-albedo feedback (e.g., Klok and Oerlemans, 2004; Thackeray and Fletcher, 2016).

The Snow-Atmosphere-Soil Transfer (SAST) model, a three-layer physical-based snow-cover forecast model, was developed in 1999 for the study of climate (Sun et al., 1999). The snow albedo in this model is parameterized with the solar angel zenith, cloud amount, snow age, and snow depth. However, the scheme does not consider the impacts of LAA on snow. In order to investigate the effects of LAA, the "Snow, Ice and Aerosol Radiation" (SNICAR) model (Flanner and Zender, 2006), a physical-based snow-albedo scheme that includes snow aging and the effects of deposited LAA, is incorporated into the SAST model. The objectives of this study are to: (2) improve the ability of the SAST model to reproduce the seasonal snow-cover process, and (3) investigate the impacts of LAA on the surface energy balance and water budget during the winter and spring months. Following this introduction, section 2 describes the development of the coupled model, including details of the SAST and SNICAR models and the aerosol layer schemes. Section 3 describes the data used for the validation and experimental design. Section 4 presents the model simulations and analyses, demonstrating the advantage of the new coupled model over the old SAST model through consideration of LAA deposition. Section 5 provides a summary and discussion.

The SAST model is a seasonal snow-cover forecast model for climate studies and hydrological use (Sun et al., 1999). This model deals with the physical processes of the mass and energy exchanges among snow, atmosphere, and soil. It contains a series of significant physical processes, such as energy balance and mass balance, which involve heat conduction, three-phase change, movement of water inside the snow, snow density, compaction, and a snowpack layer scheme. There are three main prognostic variables (enthalpy, mass, and snow depth) in the model. The snowpack is divided into three layers at most, so as to save computational resource long-term GCM integrations. A number of validations and evaluations have demonstrated that the simulation results of the SAST model are reasonable and consistent when compared with observations (e.g., Franz et al., 2008; Rutter et al., 2009; Essery et al., 2013). In addition, the model has been implemented into the Simplified Simple Biosphere Model (Xue et al., 1991) within regional and global atmospheric models (e.g., Sun and Xue, 2001; Waliser et al., 2011; Oaida et al., 2015).

According to (Gray and Landine, 1987), and (Verseghy, 1991), the snow albedo in the SAST model under clear-sky conditions is a piecewise function of snow age, which is classified by snow depth (deep or shallow) and snow temperature (dry or melting).

For melting and shallow snow cover (snow depth <25 cm), melting and deep snow cover (depth >25 cm), and for dry and cold snow cover, the snow albedo [α₀(t)] can be presented, respectively, as follows: \begin{equation} \label{eq1} \left\{ \begin{array}{rcl} \alpha_{0}(t)&=&\alpha_{0}(t-1)-\dfrac{0.071\Delta t}{86400}\\[2mm] \alpha_{0}(t)&=&0.5+[\alpha_{0}(t-1)-0.5]\exp\left(-\dfrac{0.01\Delta t}{3600}\right)\\[2mm] \alpha_{0}(t)&=&\alpha_{0}(t-1)-\dfrac{0.006\Delta t}{86400} \end{array} \right. , \ \ (1)\end{equation} where α₀(t-1) is the albedo value at the last time step, and ∆ t is the time step.

Moreover, after every new snowfall, the clear-sky albedo would increase by 0.1 for every 1 cm of new-snowfall depth, but the maximum value is fixed at 0.92.

Additionally, the snow albedo is also affected by the cloud cover and solar zenith angle, according to (Siemer, 1988). Hence, the snow albedo α(t) under cloudy conditions can be modified as shown below: \begin{equation} \label{eq2} \alpha(t)=\alpha_{0}(t)+\alpha_{0}^3(t)[1-\alpha_{0}(t)]F(N,{\rm AG}) ,\ \ (2) \end{equation} where N is the amount of cloud and AG is the minimum value between the sun elevation angle and π/3. F(N, AG) is a function of clouds amount (N) and sun elevation angle ( AG).

It can be seen that the snow-albedo process neglects the impacts of the LAA deposited on snow, and only takes semi-empirical functions to present the effects of snow depth, solar angel zenith, cloud amount, and snow age on the snow albedo.

Because snow is not opaque to the visible spectrum, after downward solar radiation is reflected by the snow surface, the absorbed solar radiation continues to penetrate through the underlying snow layers, decaying exponentially with snow depth in the snowpack. (Flanner and Zender, 2005) showed that the vertical absorption profile will not only affect the energy and mass balance of the snowpack, but also impose fluctuation in air temperature. This implies that the radiative transfer process in the snowpack plays a significant role in the simulation of snow and climate. Therefore, it is necessary to accurately simulate the radiative transfer process between snow layers.

In the SAST model, Beer's law is applied to derive the penetration of solar radiation through snow layers (Jordan, 1991), and the extinction coefficient in Beer's law depends on the snow grain diameter (d), which is parameterized as an empirical function of snow density (ρ, units: kg m^-3) as follows (Anderson, 1976): \begin{equation} \label{eq3} \left\{ \begin{array}{l@{\quad}l} d=0 & \rho>920\\[1mm] d=2.796\times 10^{-3} & 400\le \rho \le 920\\[1mm] d=1.6\times 10^{-4}+1.1\times 10^{-13}\rho^4 & \rho<400 \end{array} \right. , \ \ (3)\end{equation} As knowledge of snowpack microphysics has improved (Flanner and Zender, 2006), it has become possible to physically simulate the effects of LAA in snow, the solar radiative absorption transfer process, and snow grain size, rather than use simple empirical functions, as described in the following subsection.

The SNICAR model, developed by (Flanner and Zender, 2005) and (Flanner et al., 2007), can be applied to calculate the snow albedo and solar absorption for each snow layer within aerosols stored in the snow. The snow albedo and solar absorption vertical profiles depend on the solar zenith angle, ice effective grain size, and mass concentrations of the deposited aerosols.

Specifically, SNICAR uses a two-stream, multiple scattering, and multi-layer radiative approximation from (Toon et al., 1989), to calculate the downward and upward radiative fluxes for each snow layer, and to derive the layer absorption and the surface albedo. It requires the following bulk optical properties for each snow layer and spectral band: extinction optical depth (τ); single-scatter albedo (ω); and scattering asymmetry parameter (g). Each constituent (e.g., ice and aerosol species) in the snow layer has its own optical properties for each spectral band, which are computed offline following Mie Theory (Bohren and Huffman, 1983). Due to the complexity of the radiative transfer processes in the SNICAR model and time-consuming calculations in climate models, here, only a name list-defined lookup table is employed, which is derived offline for these optical properties according to the technical notes of CLM4.0 (Oleson et al., 2010). The lookup table provides detailed information on these optical properties in five spectral bands, including one in the visible spectrum (0.3-0.7 μm) and four in the near-infrared spectrum (0.7-1.0 μm, 1.0-1.2 μm, 1.2-1.5 μm, and 1.5-5.0 μm). The weights of each band for direct and diffuse radiation are determined with offline hyperspectral radiative transfer calculations for an atmosphere typical to midlatitude winter (Flanner et al., 2007). The optical properties contained in the lookup table obey lognormal distributions over the range of ice effective radii: 30 μm <r _e(t)<1500 μm, at a resolution of 1 μm. Hence, in SNICAR, in order to obtain these optical properties, it is imperative to firstly calculate the effective grain radius of the ice particles for each snow layer.

The evolution of the ice effective grain radius [r _e(t)] represents snow aging in the snow microphysics (e.g., Wiscombe and Warren, 1980; Flanner and Zender, 2006). r _e(t) is a function of dry-snow metamorphism (∆ r _e,dry), liquid water-induced metamorphism (∆ r _e,wet), the refreezing of liquid water (f _rfrz), and the addition of new snowfall (f _new). It can be expressed as \begin{equation} \label{eq4} r_{\rm e}(t)=[r_{\rm e}(t-1)+\Delta r_{\rm e,dry}+\Delta r_{\rm e,wet}]f_{\rm old}+r_{\rm e,0}f_{\rm new} +r_{\rm e,rfrz}f_{\rm rfrz} ,\ \ (4) \end{equation} where r _e(t-1) is the effective radius at the last time step.

∆ r _e,dry is determined by \begin{equation} \label{eq5} \Delta r_{\rm e,dry}=\left\{\left(\dfrac{dr_{\rm e}}{dt}\right)_0\left[\dfrac{\eta}{(r_{\rm e}-r_{\rm e,0})+\eta}\right]^{1/\kappa}\right\}dt , \ \ (5)\end{equation} where the parameters (dr _e/dt)₀, η and \(\kappa\) are retrieved interactively from a lookup table according to snow temperature (T), temperature gradient (dT/dz), and snow density (ρ), to save on computational resource (Oleson et al., 2010). (Flanner and Zender, 2006) showed that a warm snow temperature with a large temperature gradient and low density cause rapid snow aging, whereas snow aging slows down under a cold snow temperature, regardless of the temperature gradient and density.

∆ r _e,wet depends on the mass liquid water fraction (f _liq), as expressed below: \begin{equation} \label{eq6} \Delta r_{\rm e,wet}=\left(\frac{10^{18}C_1f_{\rm liq}^3}{4\pi r_{\rm e}^2}\right)dt ,\ \ (6) \end{equation} where C₁=4.22× 10^-13 and f _liq=w _liq/(w _liq+w _ice) represents the liquid mass fraction.

The effective grain radius of new snow and refrozen snow are fixed to r _e,0=54.5 μm and r _e,rfrz=1000 μm, respectively.

Figure 1.  Flow chart of the coupling of the SNICAR scheme into the SAST model. Inputs and outputs are presented in oval shapes. The green box indicates physical processes in the SAST model, while the orange dashed box relates to the SNICAR scheme.

The fractions f _old, f _new and f _rfrz are the ratio of old, newly fallen, and refrozen snow mass, respectively. The refrozen snow mass is defined as: if liquid water exists in the snow layer and, simultaneously, the snow temperature is below freezing point, then refreezing would happen. The value is equal to the difference in ice content between the present and previous time step.

In summary, as shown in the orange dashed box in Fig. 1, if r _e(t-1), T, dT/dz, ρ, w _liq, w _ice and f _new are given, then r _e(t) can be derived through a series of calculations. Then, the optical properties τ, ω and g of the ice particles in five spectral bands can be obtained from the lookup tables. If the snow contains LAA, the aerosol concentrations should be added as inputs as well, and then the total extinction optical length will increase its constituent part. In this study, it is assumed that the LAA particles consist of two kinds of BC (hydrophilic and hydrophobic) and four kinds of mineral dust (with diameters of 0.1-1, 1-2.5, 2.5-5, and 5-10 μm, respectively) (Oleson et al., 2010). Finally, SNICAR utilizes the two-stream radiation solution to calculate the albedo and radiative absorption at each snow layer. The results are robust and have been validated with laboratory and field measurements (Hadley and Kirchstetter, 2012).

Since the SNICAR scheme can calculate the snow albedo based on the optical properties of the ice and aerosols particles, in the coupled scheme with the SAST model, the movement of aerosols in snow layers should be accounted for.

Aerosol-particle movement processes include loading, being buried under new snowfall, exposure to air when overlying snow melts, aerosol layers emerging, and flushing with meltwater. Here are some of the assumptions that first need to be made (Oleson et al., 2010):

(1) Deposited particles are instantly mixed in the top snow layer;

(2) Inter-layer water fluxes are computed first, and then the aerosols are added;

(3) Aerosol particles are not immediately washed out before the radiative calculations are carried out;

(4) Once the mass of the LAA is flushed out with meltwater through the bottom snow layer, this part is considered as a permanent loss;

(5) When the LAA mass distribution changes along with snow-layer combination or subdivision, it is assumed that the mass concentration of the LAA is partially uniform, which means that the mass is uniformly distributed in the same snow layer, but is separate from its neighboring layers.

Changes in the aerosol mass in each snow layer due to the deposition and flushing out with melting water can be presented as follows: \begin{equation} \Delta m_i=[k(q_{i+1}c_{i+1}-q_ic_i)+D]\Delta t , \ \ (7)\end{equation} where k represents the aerosol species-dependent scavenging efficiency (Conway et al., 1996); q_i+1c_i+1-q_ic_i means the net change in flow, inflow from the overlying snow layer, and outflow to the underlying snow layer, in which q is the flow of water, c is the particle mass mixing ratio, and c_i=m_i/(w _liq,i+w _ice,i); and D is the deposition rate (only for the top layer).

In the snow model layer scheme, the top layer is very thin——typically less than 2 cm. If there is a snow storm with a new-snowfall depth of greater than 2 cm, aerosol particles will all be buried under the newly formed surface snowfall layer, and then the aerosol masses m _arsl(k,i) at the original top layer (i) are all added to the underlying snow layer (i-1), and can be expressed as follows: \begin{equation} \label{eq7} m_{\rm arsl}(k,i-1)=m_{\rm arsl}(k,i-1)+m_{\rm arsl}(k,i) .\ \ (8) \end{equation} Furthermore, the content in the newly formed top layer (i) will be determined by whether or not the new snow is clean.

During the snow ablation process, snow melting occurs. If the overlying snow layer almost melts away, but the snow depth is large enough to maintain an unchanged layer number, this indicates that the underlying snow layer will be subdivided into two layers, p and 1-p, respectively; p is the aerosol mass supplement to the top layer, and these parts of the aerosol mass of the underlying snow layer (i-1) should be added to the overlying snow layer (i), \begin{equation*} \label{eq8} m_{\rm arsl}(k,i)=m_{\rm arsl}(k,i)+m_{\rm arsl}(k,i)p , {\rm (9a)} \end{equation*} and the aerosol mass at the new sublayer is \begin{equation*} \label{eq9} m_{\rm arsl}(k,i-1)=m_{\rm arsl}(k,i-1)(1-p) . {\rm (9b)} \end{equation*} However, if the overlying snow layer (i) is not thick enough to be a separated layer, then it will merge with the underlying snow layer (i-1) into a new snow layer. Hence, all the aerosol mass at layer (i) are added to the underlying snow layer (i-1): \begin{align} \label{eq10} m_{\rm arsl}(k,i-1)&=m_{\rm arsl}(k,i-1)+m_{\rm arsl}(k,i) ;\ \ (10a)\\ \label{eq11} m_{\rm arsl}(k,i)&=0; \ \ (10b)\end{align} where m _arsl(k,i) and m _arsl(k,i-1) are the masses of aerosol species k at layer i and its underlying layer i-1, respectively.

According to (Skiles and Painter, 2017), the BC concentrations at the Swamp Angel Study Plot (SASP), Colorado, USA, during 25 March to 18 May 2013, were between four and seven orders of magnitude less than the dust concentrations. Meanwhile, utilizing laser light diffraction, (Skiles et al., 2017) took six samples, including both single and merged dust layers, to quantify the dust particle size distributions between 0.05 and 2000 μm at the SASP in spring 2013. Then, based on the measurements, the dust particles were binned into four groups (the same group classification as defined in the SNICAR model). The relative contributions were as follows: 3%, 6%, 9% and 82% for 0.1-1.0, 1.0-2.5, 2.5-5.0 and 5.0-10.0 μm, respectively. The dust particle size distributions were relatively consistent from sample to sample.

Due to a lack of records of mixed mass ratios in 2010-12, according to their (Skiles et al., 2017) measurements, we set the mixed mass ratios of BC and dust as 0.01% and 99.99%, respectively. Furthermore, the mixed mass ratios of the four kinds of dust to total dust are 3%, 6%, 9% and 82%, respectively. Still, we recognize this assignment as a limitation, because interannual variability of the dust size distributions may exist, and plan to analyze this uncertainty in the future.

Figure 1 gives an overview of the calculations performed in the coupled scheme. The physical processes in the SAST model are shown in the green box, while the orange dashed box presents the SNICAR steps. At the beginning, a series of basic physical processes in the snow, such as initialization, new snowfall, updating the snow layering, and compaction, are performed after inputting hourly meteorological data as external forcing. These processes are almost the same as in the original SAST model [further details in (Sun et al., 1999)], except the calculation of the ice grain radius in compaction is replaced later in the SNICAR model. Secondly, when the model moves on to the SNICAR part, it requires information on the effective grain radius [r _e(t-1)] from the previous time step, as well as snow temperature (T), temperature gradient (dT/dz), snow density (ρ), liquid water content (w _liq), ice content (w _ice), new-snowfall content (f _new), and the aerosol mass concentration, as inputs. After completion of the approximate radiative transfer calculations in the SNICAR model, as mentioned in section 2.2 and shown in the orange box of Fig. 1, the snow albedo and radiative absorption at each snow layer are returned to the SAST model. Then, the SAST model continues to perform the surface energy balance, and update the temperature profile, mass balance, snowmelt, water flow, snow layering, and so on. Further details concerning these processes can be found in (Sun et al., 1999). In the presence of LAA in the snowpack, the aerosol stratigraphy will also change simultaneously with the SAST snowpack layering adjustment.

The in-situ observations used in this study are from the SASP, Colorado, USA (37.9^°N, 107.7^°W, ∼3371 m in altitude). This location is highly suitable for measuring precipitation and snowpack accumulation, because the wind speeds are low enough that the redistribution of snow cover by wind is negligible. It provides the hourly forcing data, which include the air temperature, pressure, relative humidity, wind speed, precipitation, global solar radiation, downward longwave radiation, reflected radiation, and snow depth. More details regarding the meteorological data can be found in (Landry et al., 2014), and at http://snowstudies.org.

Figure 2.  Daily variations of air temperature (black line), surface temperature (blue dashed line), snow precipitation (green stars), rain precipitation (black dots), and aerosol loading mass (orange bars) at the SASP in 2010-12.

Once aerosols are deposited on snow, they will be recorded as a "dust-on-snow" event (the dominant aerosol at the SASP is dust). For each dust-on-snow event, the dust loading is determined by collecting the dust layer and some clean snow above and below the dust layer in a column over a 0.5-m² area. After the samples are melted, dried, and weighed, the dust loading mass flux of the event is recorded (g m^-2) (Painter et al., 2012). The dust loading time and mass flux at the SASP in 2010-12 are presented in Fig. 2 (orange bars).

For each dust-on-snow event, after it is blanketed by the new snowfall later on, a discernible mixed dust-in-snow layer is formed, and the dust remains in the layer in which it was originally deposited (e.g., Conway et al., 1996; Flanner et al., 2007). The position of each dust-on-snow event is recorded in regular snow-pit measurements. Combined with the measurements of each dust loading mass and its corresponding position, the vertical profile of the dust mass flux in the snow can be derived. Measurements from snow pits at the SASP are sampled monthly during the snow accumulation period, and weekly during the snow ablation period. Dust mass fluxes in the top layer are estimated from the vertical profiles to validate the coupled model's aerosol simulations. More details regarding the aerosol data can be found in (Painter et al., 2012), and at http://www.codos.org/#codos.

Figure 2 displays the time series of air temperature (black line), surface temperature (blue dashed line), snow precipitation (green stars), rain precipitation (black dots), and aerosol loading mass (orange bars) at the SASP in 2010-12. It can be seen that the air temperature is very low, with the minimum value reaching -26.9^°C and -17.9^°C in 2010/11 and 2011/12, respectively. Snow cover lasts seven to eight months per year at the SASP, and the daily maximum total snow precipitation reaches 54.0 kg m^-2 and 46.4 kg m^-2 in 2010/11 and 2011/12, respectively. Interestingly, there are interannual variations in spring temperature and precipitation between 2011 and 2012. From March to May, the average air temperature is -1.9^°C in 2011. In comparison, the value is 0.7^°C in 2012, which is much higher than that in 2011. In addition, from March to May in 2011, persistent snowfall hits the SASP, and the seasonal total precipitation is as high as 516 kg m^-2, whereas the value is only 153 kg m^-2 in 2012. These differences produce dramatically opposing snowmelt scenarios between 2011 and 2012. Specifically, in 2011, larger amounts of new snowfall and lower temperatures produce a very large snowpack and very late dates of snow ablation. By contrast, in 2012, less snowfall and sunny, sometimes unseasonably warm weather, begins in early March, enabling snowpack ablation to begin early. However, aerosol loadings are almost all occur during late-winter and spring, with the maximum loading happening on 21 March in 2011 (5.124 g m^-2) and on 18 March in 2012 (3.277 g m^-2).

The experiments are conducted under three scenarios: (2) SAST _original (original SAST scheme without any change in the SAST); (3) SAST _pure (SNICAR scheme coupled with the SAST model but the snow is clean); (4) SAST _arsl (coupled model and dirty snow; aerosol loading data used as the forcing data). By comparing Scenario 1 and Scenario 3 against observations, we can validate the coupled model. In addition, by comparing Scenario 2 and Scenario 3, we can evaluate the impacts of the LAA on the snowpack energy balance and water budget.

In this section, meteorological forcing data and the aerosol loading data at the SASP during 2010-12 are used to validate the SAST _arsl.

According to the SNICAR scheme, it is the mass concentrations and optical properties of ice particles and aerosol particles that mainly determine the snow albedo and radiative absorption at each snow layer. Therefore, a proper simulation of aerosol mass is necessary. Before comparing the aerosol mass in the top snow layer, it is important to keep in mind that the snow surface is not fixed. The aerosol mass content in the top snow layer will increase instantly when aerosols are deposited on the snow surface. During early spring, there is still new snowfall at the SASP intermittently. The aerosol layer surface will usually be buried in the underlying snow layer by the following new snowfall. If the new snowfall is clean, then the aerosol content in the top snow layer will plunge to zero, and the content will increase again if another deposition happens——a pattern that will repeat until no new snowfall follows. Finally, during snowmelt, the aerosol content in the top snow layer will gradually increase as the underlying aerosol layers emerge.

Figure 3 compares the aerosol mass simulations between the SAST _arsl in the top snow layer and the observations in the top 3 cm of the snowpack. The dominant aerosol at the SASP is dust; hence, in this section, this is the aerosol type to which we mainly refer. It can be seen that the SAST _arsl simulations represent the changing dynamics of the aerosol mass flux in the top snow layer, with many peaks in both 2011 and 2012. Yet, the observations are quite flat in early-spring, as constrained by the temporal resolution (weekly) of the observation. During early ablation, the SAST _arsl simulation underestimates the aerosol mass, whereas during late ablation the simulation overestimates slightly. As shown in Table 1, in 2011 (2012), the average bias is -0.349 g m^-2 (0.453 g m^-2), the root-mean-square error (RMSE) is 0.946 g m^-2 (0.688 g m^-2), and the correlation coefficient is 0.982 (0.823). The above results suggest that the aerosol stratigraphy scheme is effective and reliable, thus laying a good foundation for the simulation of the snow albedo and snow depth.

Figure 4 compares the snow albedo between SAST _original and SAST _arsl with respect to the observations. The results show that, in general, the variation of the snow albedo in SAST _original fits well with the observations during the snow accumulation period, but fails to reproduce the snow albedo value during snow ablation. This may be related to the fact that SAST _original does not consider the radiative forcing of LAA. Previous studies have demonstrated that, apart from snow aging, aerosol radiative forcing is also a significant contributor to albedo variations (Warren and Wiscombe, 1980). The average bias for the snow albedo in the SAST _original simulation is 0.097, as shown in Table 1.

Figure 3.  Comparison of aerosol mass between simulations by SAST _arsl (blue line) in the top snow layer and observations in the top 3 cm (red diamonds). Brown reference lines are marked as the dates when aerosol loading occurs. The simulations are the results of the daily mean of the model outputs (hourly).

Figure 4.  Comparison of snow albedo between simulations by SAST _original (green line) and SAST _arsl (blue dashed line) against observations (red line) at the SASP in 2011 and 2012. Brown bars are marked as the dates when aerosol loading occurs.

In contrast, irrespective of whether in the snow accumulation or snow ablation period, the values in SAST _arsl (blue dashed line) are closer to the observations than those in SAST _original, except for the maximum values. The simulation of snow albedo by SAST _arsl is closely dependent on the distribution of aerosol mass. As mentioned above, the aerosol mass simulation is slightly smaller than the observations during the later snow ablation period, so the modeled albedo is relatively higher than observed. As displayed in Table 1, the average bias for snow albedo in 2011 (2012) according to SAST _arsl is -0.006 (0.024), the RMSE is 0.061 (0.099), and the correlation coefficient is 0.917 (0.889). The average RMSE for the SAST _arsl has decreased by around 6%.

The simulations of snow depth in SAST _arsl match well with the observations. Figure 5 compares the snow depth simulated in SAST _original and SAST _arsl against the observations. The results are consistent with the aerosol mass and snow albedo. Due to the exclusion of the aerosol radiative effect on snow cover, SAST _original overestimates the snow depth at the SASP during the snow ablation period, with an average bias of 11.4 cm higher than observed. In contrast, SAST _arsl represents the evolution of snow cover better than SAST _original. In 2011 and 2012, the average bias is 0.5 cm and -1.5 cm, respectively, and the RMSE is 13.7 cm and 8.3 cm, which demonstrates an improvement compared with SAST _original.

From the above analyses, we can conclude that, to obtain a better simulation, the aerosol radiative transfer process should be included in snow albedo scheme in the SAST model. SAST _arsl can be used reasonably and reliably in climate studies to obtain a proper simulation of the aerosol mass distribution, snow albedo, and snow depth.

Figure 5.  Comparison of snow depth between simulations by SAST _original (green line) and SAST _arsl (blue dashed line) against observations (red line) at the SASP in 2011 and 2012. Brown bars are marked as the dates when aerosol loading occurs.

In this section, the impacts of aerosols are further investigated by comparing the simulation results between clean and dirty snow.

4.2.1. Albedo and snow depth

From Fig. 6, it can be seen that, during the snow accumulation period, there is little difference in the simulated snow albedo between SAST _pure and SAST _arsl. The snow albedo either increases due to the new snowfall, or decreases due to snow aging. The impact of the LAA cannot be investigated in this period, because the aerosol mass content in the top layer is usually quite low at the SASP (Fig. 3). Most of the aerosol loading happens in late-winter and spring, and some already-deposited aerosol particles are usually being buried in the underlying snow layer by the following new snowfall. In the snow ablation period, the results of SAST _pure are markedly different from the observations. In contrast, SAST _arsl can capture the features of the albedo's variations, and the results match well with the observations, as error analyses described in Table 1. Similar conclusions are also drawn by comparing the simulations of snow depth in the two scenarios (Fig. 7). Therefore, the comparisons of the snow albedo and snow depth simulations between SAST _pure and SAST _arsl demonstrate the significance of including the effects of the deposited LAA in snow, especially during snow melting. Snow aging alone cannot fully explain the decrease in snow albedo during snow ablation.

4.2.2. Aerosol effects on the energy balance

Figures 8c and d compare the effective grain size at the top snow layer between SAST _pure and SAST _arsl. In general,during the snow accumulation period, the snow effective grain size remains at a small value, below 150 μm. By contrast, during the snow ablation period, the snow effective grain size increases gradually. This is partially due to snow metamorphism. On the other hand, as shown in Figs. 8e and f, the effective grain size in SAST _arsl is greater than in SAST _pure, and the time corresponds well with aerosol existing in the top snow layer, as shown in Figs. 8a and b. This implies that the presence of aerosol particles is favorable for the growth of snow grain size, which further enhances radiative absorption; the larger the snow grain size, the more radiative absorption there is by ice particles. This additional radiative absorption is known as the "aerosol first indirect radiative forcing effect" (Hansen and Nazarenko, 2004).

Figure 6.  Comparison of snow albedo between simulations and observations (red line) at the SASP in 2011 and 2012. The blue dashed line is SAST _arsl, meaning snow containing aerosols. The green line is SAST _pure, meaning clean snow. Brown bars are marked as the dates when aerosol loading occurs.

Figure 7.  Comparison of snow depth between simulations and observations (red line) at the SASP in 2011 and 2012. The blue dashed line is the simulation by SAST _arsl, meaning snow containing aerosols. The green line is SAST _pure, meaning clean snow. Brown bars are marked as the dates when aerosol loading occurs.

Figure 8.  Comparison between SAST _pure (purple line) and SAST _arsl (blue dashed line) in terms of (a, b) the aerosol mass distribution at the top snow layer (blue line) against observations (red diamonds), where brown bars are marked as the dates when aerosol loading occurs; (c, d) the effective radius at the top snow layer; (g, h) the radiative absorption at the top snow layer; and (i, j) runoff. Panels (e, f) show the difference in effective radius between SAST _pure and SAST _arsl at the top snow layer (blue line). Panels (a, c, e, g, i) are the results for 2011, while (b, d, f, h, j) are for 2012.

Figures 8g and h compare the radiative absorption at the top snow layer between SAST _pure and SAST _arsl. The snow cover containing aerosols (blue dashed line) absorbs more solar radiative energy than the clean snow (purple line), especially when the aerosol concentration in the snow layer is high. Considering aerosol loading as an external forcing factor makes the aerosol itself increase radiative absorption directly, its first indirect radiative forcing effect cannot be neglected. Moreover, the extra energy accelerates the melting rate of the snow cover, which leads to exposure of the dark substrate being much earlier than usual, and thus brings more energy to the earth system. This is thought to be the second indirect radiative forcing effect of aerosol (Hansen and Nazarenko, 2004).

To quantify the impacts of the LAA on the surface energy balance, the radiative forcing is calculated. As shown in Fig. 7, we mark the time as Phase 1, which is defined as the period from 20 March to the "snow-all-gone day" in SAST _arsl. Next, Phase 2 is defined as the period from the snow-all-gone day in SAST _arsl to the snow-all-gone day in SAST _pure. The radiative forcing during Phase 1 is the sum of the aerosol direct effect and first indirect effect; while during Phase 2, the radiative forcing is equal to the aerosol second indirect radiative forcing effect only. As shown in Table 2, during Phase 1, the radiative forcing is 20.6 W m^-2 and 26.1 W m^-2 in 2011 and 2012, respectively, and 23.3 W m^-2 on average. During Phase 2, the values are much larger, at 132.7 W m^-2 in 2011 and 146.4 W m^-2 in 2012, and 139.6 W m^-2 on average. During the whole ablation period, the radiative forcing is 39.9 W m^-2 in 2011 and 55.0 W m^-2 in 2012, and 47.5 W m^-2 on average. (Painter et al., 2007) and (Skiles et al., 2012), using semi-empirical methods to calculate the second indirect radiative forcing at the SASP in 2006, reported values of 147 8 W m^-2 and +150 W m^-2, respectively; our results are consistent with their studies.

4.2.3. Aerosol effects on the water budget

Extra radiative absorption due to the existence of LAA in snow accelerates the snow melting rate and advances the timing of the snow-all-gone day, which alters the amount and timing of streamflow. Table 3 lists the snow-all-gone days for SAST _pure and SAST _arsl with respect to observations. It can be seen that, on average, the LAA cause the snow-all-gone day to advance 19.5 days. At the same time, the peak runoff discharge and timing advance by 18 days, compared to the clean-snow scenario, which is consistent with the snow-all-gone day. These changes will have substantial implications for the hydrologic cycle and atmospheric processes.

Snow albedo is well-known as a crucial parameter of snow for modulating the energy exchange among snow, soil, and atmosphere. With the increase in aerosol emissions due to human activity, the effects of LAA in snow cannot be ignored in snow and climate studies. In this paper, the original snow albedo scheme and radiative transfer scheme in the SAST model are replaced by the SNICAR scheme, which is a physical-based albedo scheme that includes snow aging and the effects of deposited LAA in snow. Furthermore, aerosol-particle movement, including loading on the snow surface, being buried under new snowfall, exposure to the air when the overlying snow melts, aerosol layers emerging, and flushing with meltwater, are also described in the model.

Compared to the empirical-function calculations in the original SAST, the coupled model presents a major step in model development towards improving snow-albedo simulation based on the physical-based radiative transfer process and aerosol-in-snow effects. In summary, the model can be used to:

(1) Physically and realistically simulate the snow albedo;

(2) Reproduce the evolution of seasonal snow cover, especially during the snow ablation period;

(3) Present the evolution of effective snow grain size, allowing for snow metamorphism processes and the presence of aerosols in the snow;

(4) Simulate aerosol-particle movement processes and mass distribution along with snow-layer adjustment;

(5) Evaluate the impacts of LAA on the surface energy balance and water budget.

Measurements at the SASP during 2010-12 are used to validate the coupled model. The results show that the depiction of aerosol movement is quite good, as the modeled aerosol mass content in the top snow layer is generally consistent with observations. The SNICAR scheme, which contains snow metamorphism and the radiative effects of deposited LAA in snow, enables SAST _arsl to physically mimic the evolution of snow grain size. As a result, simulations of the variation in snow albedo and snow depth by SAST _arsl are better than those of SAST _original, especially during the snow ablation period. This implies that the coupled model, SAST _arsl, can be used to properly simulate aerosol mass and the evolution of snow albedo and snow depth.

The impacts of LAA on the snowpack are also investigated, through comparisons between simulations of clean snow and dirty snow. Radiative forcing and runoff are evaluated to illustrate the impacts of LAA on the surface energy balance and water budget, separately. On average, the total radiative forcing of aerosol loading on the seasonal snow cover at the SASP is 47.5 W m^-2. The snow-all-gone day advances by 19.5 days, thus altering the timing and amount of peak runoff. These changes will have a profound influence on the following hydrological cycle and atmospheric processes.

As emissions of aerosols are increasing in the context of climate change, snow cover is increasingly at risk of LAA deposition. To better present snow cover in land-surface and climate models, apart from the meteorological and solar radiation data, it is also crucial to include aerosol loading data as an external forcing input, and incorporate aerosol radiative transfer processes into the snow-albedo scheme. Additional measurements and further analyses are needed for model development and validation.

It is important to note that several additional improvements are still needed. For example, it would be beneficial to extend SAST _arsl for regional climate studies. Besides, when aerosols emerge and enhance the melting process, the surface roughness will increase at the same time. This additional aerosol-induced surface roughness may enhance the sublimation rate of ice and the evaporation rate of meltwater. Unfortunately, these processes and their effects on snowmelt runoff are still unclear. Measurements and modeling are needed to better understand these effects.