The role of surface energy fluxes in pan-Arctic snow cover changes

We analyze snow cover extent (SCE) trends in the National Oceanic and Atmospheric Administration’s (NOAA) northern hemisphere weekly satellite SCE data using the Mann–Kendall trend test and find that North American and Eurasian snow cover in the pan-Arctic have declined significantly in spring and summer over the period of satellite record beginning in the early 1970s. These trends are reproduced, both in trend direction and statistical significance, in reconstructions using the variable infiltration capacity (VIC) hydrological model. We find that spring and summer surface radiative and turbulent fluxes generated in VIC have strong correlations with satellite observations of SCE. We identify the role of surface energy fluxes and determine which is most responsible for the observed spring and summer SCE recession. We find that positive trends in surface net radiation (SNR) accompany most of the SCE trends, whereas modeled latent heat (LH) and sensible heat (SH) trends associated with warming on SCE mostly cancel each other, except for North America in spring, and to a lesser extent for Eurasia in summer. In spring over North America and summer in Eurasia, the SH contribution to the observed snow cover trends is substantial. The results indicate that ΔSNR is the primary energy source and ΔSH plays a secondary role in changes of SCE. Compared with ΔSNR and ΔSH, ΔLH has a minor influence on pan-Arctic snow cover changes.


Introduction
The pan-Arctic domain is one of the most sensitive regions on Earth to global climate change (Manabe and Stouffer 1994, Miller and Russell 2000, Holland and Bitz 2003, Serreze et al 2009. As the largest single component of the cryosphere in terms of spatial extent (Armstrong and Brodzik 2001), snow cover variations in space and time have resulted in significant changes in the surface energy and water budgets over the pan-Arctic land region (Serreze et al 2000, Peterson et al 2002  the basin scale in two regions of Alaska (Robinson 1986); the Red River Valley of North Dakota and Minnesota (Dyer and Mote 2002); Trail Valley Creek of northern Canada Pomeroy 1996, Pohl andMarsh 2006). However, few studies have examined the large-scale factors that would provide better understanding of the relative importance of snow surface energy balance components (Male and Granger 1981, Cline 1997, Leathers et al 2004, in large-scale snow cover changes, especially for the pan-Arctic land area. On the other hand, land surface models (e.g. Liang et al 1994) have improved to the point that they may, in some cases, serve as surrogates for in situ observations. Off-line runs of these models provide sources for most or all terms in the snow surface energy balance and offer the opportunity to investigate the nature of the space-time variability of the snow surface energy budget (Betts et al 2009, Troy and Wood 2009, Shi et al 2010. In this letter, our main objective is to identify the individual role of surface energy fluxes in the snow surface energy balance and determine which are most responsible for observed changes in snow cover extent (SCE) of the pan-Arctic land region as shown in figure 1. First, monotonic trends in satellite snow cover observations and corresponding reconstructions generated by a land surface model are analyzed using a non-parametric trend test. Subsequently, the relationships between observed SCE and modeled surface radiative and turbulent fluxes are examined. In the final part of the letter, the relative importance of each component in the snow surface energy balance is estimated.

Snow cover extent data
Observed seasonal values of SCE were extracted from the weekly snow cover and sea ice extent version 3 product for the northern hemisphere maintained at the National Snow and Ice Data Center (NSIDC), which combines snow cover and sea ice extent for the period from October 1966 through June 2007 (Armstrong and Brodzik 2005). The data set is based on weekly maps of continental SCE produced by the National Oceanic and Atmospheric Administration's (NOAA) National Environmental Satellite Data and Information Service (NESDIS) (Robinson et al 1993, Frei andRobinson 1999), which were derived from digitized versions of manual interpretations of Advanced Very High Resolution Radiometer (AVHRR), Geostationary Operational Environmental Satellite (GOES), and other visible band satellite data. This satellitebased data set has been regridded to the NSIDC EASE grid with a spatial resolution of 25 km. Our study is restricted to the period 1972-2006 because there are some missing charts between 1967 and 1971 (Robinson 2000). Although ending the time series in 2006 leaves out some exceptionally low Arctic spring SCE values in recent years (e.g. 2008-10), the nonparametric statistical method we used (section 3.1) is robust to modest changes in the length of the record analyzed. In addition, Greenland is not included in the analyses as its snow cover is mainly perennial in nature (Déry and Brown 2007). Brown et al (2010) have assessed this SCE record (commonly referred to as the NOAA weekly SCE record) in comparison to other available Arctic snow cover data sets. In general, their study, and others (Wiesnet et al 1987, Robinson et al 1993 have found that the NOAA weekly SCE data set is reliable for continental-scale studies of snow cover variability. It has become a widely used tool for deriving trends in climate-related studies (Groisman et al 1994, Déry and Brown 2007, Flanner et al 2009, Derksen and Brown 2011, notwithstanding uncertainties in some parts of the domain for certain times of the year, especially during springtime over northern Canada (Wang et al 2005). A more recent update to the data set we used (NOAA snow chart climate data record (CDR)) is now available (Brown and Robinson 2011), but the differences between the new CDR and the data set we used at the pan-Arctic scale are small.
In this study, simulated SCE was reconstructed from 1972 to 2006 using the variable infiltration capacity (VIC) model, which is a macroscale hydrologic model that solves the energy and water balance and represents ephemeral snow cover over a gridded domain (Liang et al 1994(Liang et al , 1996. The off-line simulations from VIC used here are at a 3 h time step in full energy balance mode (meaning that the model closes its surface energy budget), forced with daily precipitation, maximum and minimum temperatures and wind speed through 2007 at a spatial resolution (EASE grid) of 100 km, constructed using methods outlined by Adam and Lettenmaier (2008). Precipitation and temperature were from gridded observations (Willmott and Matsuura 2009) and wind speed was from the NCEP/NCAR reanalysis (Kalnay et al 1996). Precipitation was adjusted for gauge undercatch and orographic effects as described by Adam and Lettenmaier (2003) and Adam et al (2006). The snow parameterization in VIC represents snow accumulation and ablation processes using a two-layer energy and mass balance approach (Cherkauer and Lettenmaier 2003) and a canopy snow interception algorithm (Storck et al 2002, Andreadis et al 2009 when an overstory is present. In the VIC model, each grid cell is partitioned into five elevation (snow) bands, which can include multiple land cover types (tiles). The snow model is then applied to each tile separately. When snow water equivalent is greater than zero, VIC assumes that snow fully covers the tile. For each grid cell, the simulated snow cover extent is calculated as the area averages of the tiles.

Surface energy fluxes data
Surface energy fluxes including downward shortwave radiation (DSW) and downward longwave radiation (DLW) were calculated by using a Temperature INDex (TIND) scheme (Kimball et al 1997, Thornton and Running 1999, Shi et al 2010 wherein DSW and DLW are estimated based on relationships with the daily temperature range and daily average temperature, respectively. TIND has been commonly used in model intercomparison experiments such as the Project for Intercomparison of Land Parameterization Schemes (PILPS) (e.g. Pitman et al 1999) and land surface models, such as VIC, for long-term simulations in cases when direct observations of energy fluxes are not available. Shi et al (2010) evaluated DSW, DLW, and albedo computed in an off-line simulation of VIC embedded with TIND along with satellite data and global reanalysis products in comparison with in situ observations from the Global Energy Balance Archive (GEBA, Ohmura et al 1989) and showed that these estimates compared well with observations over the pan-Arctic land region. Compared to the in situ observations, the mean seasonal DSW from the European Centre for Medium-Range Weather Forecast (ECMWF) 40 Year Reanalysis (ERA-40), the ECMWF Interim Reanalysis (ERA-Interim), the International Satellite Cloud Climatology Project Flux Data (ISCCP-FD), and the TIND-based scheme all have small biases (±20 W m −2 ). ERA-40, ERA-Interim and VIC DLW deseasonalized monthly anomalies had high correlations (r = 0.96, 0.97 and 0.91, respectively) with GEBA observations whereas the correlation for the satellite-based (ISCCP-FD) product was somewhat lower. VIC deseasonalized monthly albedo had similar anomaly correlations with GEBA observations as did ERA-40, ERA-Interim and ISCCP-FD estimates.
Surface net radiation (SNR) was obtained as the sum of net shortwave (SW) and longwave radiative (LW) fluxes. SW at the snow surface is a measure of the difference between DSW and upward shortwave radiation (USW). USW is the product of DSW and snow surface albedo that is assumed to decay with age based on relationships published by the US Army Corps of Engineers (1956). LW is the sum of DLW emitted by the atmosphere and the fluxes emitted upward by a melting snow surface. DLW was estimated using equation (2.42) from Bras (1990), which is based on air temperature and a function for emissivity from Tennessee Valley Authority (1972). The turbulent fluxes (sensible heat (SH) and latent heat (LH)) near the snow surface were produced using VIC's bulk aerodynamic approach, which is described in Andreadis et al (2009). The algorithm requires snow surface temperature, which is calculated by VIC's snow algorithm, wind speed, surface air temperature and relative humidity, the last three of which are taken from the forcing data. Bulk transfer coefficients for momentum, heat, and water vapor were calculated initially for a neutral condition (Price and Dunne 1976). Subsequently, the aerodynamic resistance in the presence of snow cover was corrected using the bulk Richardson's number for stable or unstable atmospheric conditions (Anderson 1976) as implemented in the VIC snow model (Andreadis et al 2009). A similar approach has been successfully applied in various settings in Arctic environments (e.g. Hinzman et al 1991, Woo et al 1999, Boike et al 2003. Other energy fluxes in the snow surface energy balance, such as ground heat and advective fluxes were also generated by VIC. VIC's energy flux convention is that surface energy fluxes toward the snow surface are defined as positive.

Temporal analyses of SCE and surface energy fluxes
We defined April and May as spring, consistent with Groisman et al (1994), and a three-month window centered on July was defined as summer. Observed and simulated SCE in spring and summer expressed as a snow-covered fraction were calculated for the Eurasia and North America study domains (figure 1). Figure 2 shows North American and Eurasian snow cover fraction (SCF) during spring and summer from VIC and NOAA observations for the period of 1972-2006. The VIC simulations match the observed SCF over both North America and Eurasia quite well, with a mean absolute bias of 4.5% in spring and 0.6% for summer.
To examine long-term trends in SCE and surface energy fluxes, we used the non-parametric Mann-Kendall trend test (Mann 1945) for trend significance, and the Sen method (Sen 1968) to estimate trend slope. A 5% significance level (two-sided test) was specified. Trend tests were performed for seasonal SCE from VIC and satellite observations, and modeled radiative and turbulent fluxes averaged over the snowcovered portions of Eurasia and North America, respectively. Table 1 summarizes trend test results for spring and summer SCE over North America and Eurasia from VIC and satellite observations. Strong negative trends were found in the NOAA satellite observations, during spring and summer in both North American and Eurasian sectors of the Arctic, which are statistically significant ( p < 0.025) except that for Eurasia in spring the significance level is p < 0.10. VIC reproduces the same trend directions, with similar significance levels as compared with the satellite observations, for both continents in spring and summer. However, the VIC and observed trend slope magnitudes differ somewhat, especially in summer. The simulated trend slope in summer is about one third of that observed in North America, and 40% of that observed in Eurasia in summer. This discrepancy may be related to uncertainties in the NOAA SCE data set in July and August for both continents (Déry and Brown 2007) and in May and June for North America (Wang et al 2005). Because the trend directions and statistical significances are consistent, we chose to include July and August to maintain completeness of the summer period.
The non-parametric Mann-Kendall trend test was also applied to the surface energy inputs to the snow surface, including SNR, SH and LH. Table 2 summarizes their trends over both continents during spring and summer. Trends were computed for seasonal mean fluxes over the snow-covered portion of the regions. Strong positive trends were found in SNR during spring and summer in North America and Eurasia. Similarly, SH fluxes also had statistically significant upward trends, except for Eurasia in spring for which the direction of change was positive, but the trend was not significant at p < 0.10. LH changes were mostly negative in spring to summer, but were statistically significant at p < 0.10 only for North America in summer and Eurasia in spring.

Correlations between observed SCE and modeled surface energy fluxes
The Pearson's product moment correlation coefficient was used to assess relationships between satellite SCE and VICsimulated surface energy fluxes. Figure 3 shows scatterplots of the surface energy fluxes (a) SNR, (b) SH, and (c) LH from VIC against observed SCF over North America and Eurasia in spring and summer. The correlations of surface energy fluxes with SCF as shown in figure 3 are all statistically significant at p < 0.025 (two-sided test). The negative sign indicates that SNR and SH have opposite change directions with SCF. SNR has a relatively stronger relationship with SCF than SH, especially in Eurasia. For LH, the correlations are positive and are stronger in North America than in Eurasia.

Role of surface energy fluxes in snow cover changes
In general, the energy balance at a snow surface includes net radiative fluxes, sensible and latent heat fluxes, ground heat fluxes, and the energy transfer due to rain on snow. Over the pan-Arctic, the ground heat flux is a small component of the energy balance of a melting snowpack compared with radiative and turbulent heat fluxes. Therefore, its effects on total snowmelt can safely be ignored (Gray and Prowse 1993). Similarly, rain on snow has an important influence on the water retention characteristics of snow and water movement in the pack but is of minor importance compared with other energy fluxes (Male and Granger 1981). Surface energy fluxes toward the snow surface are defined as positive, therefore the net radiative and sensible heat fluxes usually have positive sign and supply the energy available for snowmelt. Latent heat fluxes are directed away from the snow surface and reduce the melt energy. However, it was not clear which individual component(s) of the snow surface energy budget dominates the significant spring and summer SCE recession since 1972 as shown in table 1. Therefore, we calculated the increment due to monotonic trends in SNR, SH and LH over the 35 yr period as SNR, SH and LH in order to determine the role of each term in the observed downward trends of SCE.  Table 3. Trend analyses for the terms related with SNR during spring and summer over snow-covered North America and Eurasia generated from VIC. The significance level ( p-value) was calculated by two-sided Mann-Kendall trend test. (The unit of total change during the 35 yr period for fluxes is W m −2 ; the unit for T min and T max is • C; no unit for Albedo and RH.)  Figure 4 shows the relative role of surface energy fluxes averaged over the snow-covered portions of North America and Eurasia in spring and summer. It is apparent that SNR is the dominant energy source in both spring and summer, accounting for between 61.8% and 102.3% of the energy available for snow cover changes. The contribution of SH plays a secondary role (from 25.3% to 50.0%). LH is always opposite in sign with SH and SNR and almost completely cancels SH in North America during summer and Eurasia in spring. However, LH has a smaller absolute value than SH at other times, such as in North America during spring (−11.8%) and Eurasia in summer (−10.8%). Therefore, we conclude that LH has a minor influence on pan-Arctic snow cover changes compared with SNR and SH. Figure 5 summarizes the latitudinal variations in spring and summer surface energy fluxes over North America and Eurasia. Basically, SNR has a latitudinal pattern, which in general is at a maximum in the lower latitudinal band, and then decreases with latitude poleward. Compared with SNR, LH shows a similar pattern with a negative sign in most cases. However, it sharply decreases to zero at higher latitudes. In Eurasia during summer, it even has a positive sign for the bands 70-75 • N and 75-80 • N.
Corresponding to the patterns in SNR and LH, SH is variable depending on the seasons and latitudinal bands. As shown in figure 5(a), the contribution of SH is much larger than SNR for the 45-50 • N latitudinal band in North America, reaching 78% of the total energy attributable to SCE changes. Thereafter, this contribution decreases gradually with latitude to the 70-75

Summary and discussion
By exploring long-term trends in satellite observations of SCE, we have shown that North American and Eurasian snow cover over the pan-Arctic declined significantly in spring and summer for the period 1972-2006. Furthermore, longterm means of seasonal SCE and their trend directions are reproduced by the VIC model, which allowed us to diagnose the causes of the observed trends. We have also shown that surface radiative and turbulent heat fluxes simulated in VIC have strong correlations with observed SCE. We find that positive trends in SNR are mostly associated with the observed and model-derived SCE trends. Modeled LH and SH trends associated with warming mostly cancel, except for North America in spring, and to a lesser extent for Eurasia in summer, when the SH contribution to the SCE trends remains substantial. Our results indicate that SNR is the primary energy source and SH plays a secondary role in changes of SCE. Compared with SNR and SH, LH only has a minor influence on pan-Arctic snow cover changes.
Changes in sensible and latent heat fluxes are mostly dominated by increases in pan-Arctic surface air temperature, which have risen at a rate almost twice as large as the global average in recent decades (Lugina et al 2006, Serreze and Francis 2006, Bekryaev et al 2010. As shown in table 3, the increases in SNR are mainly associated with increased SW and increased DLW due to warmer atmospheric temperature, whereas emitted upward longwave fluxes do not change much as the snowpack temperature is mostly isothermal during the melt period. Strong upward trends in SW mostly result from statistically significant decreasing trends in snow surface albedo, while the contribution from increased DSW trends is minor (<5%) except for North America during spring when  DSW decreases. VIC was forced by DSW calculated using the method of Thornton and Running (1999) based on the daily temperature range and vapor pressure. Therefore, there is a decreasing DSW trend in table 3 for North America during spring because daily minimum temperature (T min ) has increased more rapidly than daily maximum temperature (T max ) and the relative humidity (RH) has decreased in the mean time.
Aside from the reasons described above, the difference in surface energy fluxes between North America and Eurasia may be related to the data used to force VIC, which generally are of higher quality in North America than Eurasia because of denser observation networks (Niu and Yang 2007), especially for the snow-covered portions during summer. Therefore, further study is warranted to investigate the effects of variations in forcing variables such as precipitation, wind speed and cloud cover (not used in the Thornton and Running method), which affect both surface radiative and turbulent heat fluxes. Furthermore, our study has investigated only aggregate changes over North America and Eurasia; further study is warranted to resolve spatial variations in trends over these large subcontinental areas.