Spring snowmelt onset has occurred earlier across much of the Northern Hemisphere land area in the last four decades. Understanding the mechanisms driving spring melt has remained a challenge, particularly in its spatial and temporal variability. Here, melt onset dates (MOD) obtained from passive microwave satellite data are used, as well as energy balance and meteorological fields from NASA’s Modern-Era Retrospective Analysis for Research and Applications, to assess trends in the MOD and attribute melt onset across much of Arctic and sub-Arctic Eurasia and North America during the spring snowmelt season from 1979 to 2012. Across much of the Northern Hemisphere MOD has occurred 1–2 weeks earlier over this period, with the strongest trends in western and central Russia and insignificant trends across most of North America. Trends in MOD are reflected by those in energy balance terms, with energy advection providing an increasing proportion of melt energy in regions with the strongest MOD trends. Energy advection plays a larger role in melt onset in regions where snow begins melting in March and April, while insolation and longwave radiation drives melt where the MOD occurs in May and June. This implies that there is a potential shift in snowmelt drivers toward those involved in advective processes rather than radiative processes with an earlier MOD. As the high latitudes warm and terrestrial snow cover continues to melt and disappear earlier in the spring, it is valuable to elucidate regional snowmelt sensitivities to better understand regional responses to changing climatological processes.
The Arctic climate system has undergone rapid changes in the late twentieth and twenty-first centuries (Lemke et al. 2007). Increased global temperatures are associated with accelerating losses in sea ice, glacial mass balance, permafrost, and snow cover extent and duration (Serreze et al. 2009; Liston and Hiemstra 2011; Camill 2005; Comiso et al. 2008; Johannessen et al. 2004; Post et al. 2009). Amplified warming in the Arctic relative to lower latitudes has occurred primarily as a result of feedbacks in the cryosphere, particularly the reduction in albedo associated with spring snow loss (Groisman et al. 1994; Serreze and Francis 2006; Déry and Brown 2007). This amplified warming also potentially influences midlatitude weather patterns (Francis and Vavrus 2012). Snow cover is a major component of the cryosphere and earth’s climate system, playing a central role in some of these feedbacks. In addition, the spring snowmelt season is an important element of hydrologic and ecological systems (Nijssen et al. 2001; Barnett et al. 2005) and also serves as a measure of change across the high latitudes. However, despite the importance of terrestrial snow cover, spatiotemporal variability in spring snowmelt drivers is poorly understood.
As the Arctic warms, terrestrial spring snowmelt (particularly snow disappearance) has occurred an average of two to four weeks earlier than it did in the late 1970s (Foster et al. 2008; Tedesco et al. 2009; Takala et al. 2009; Wang et al. 2013), although attributing these changes has been a challenge (Rupp et al. 2013). June snow cover extent has declined the most rapidly during this period, faster than the trend in September sea ice loss (Derksen and Brown 2012) and by nearly 50% since 1967 (Brown et al. 2010). Northern Hemisphere snow cover across the entire spring season has shown similar earlier melt trends responding to warmer temperatures and changes in atmospheric circulation (Dye 2002; Brown 2000; Déry and Brown 2007), with less pronounced trends over North America than Eurasia (Dyer and Mote 2006; Brown and Robinson 2011). The regional and monthly differences in these trends suggest that melt drivers may exhibit considerable variability, requiring attribution that adequately resolves these differences.
It is well understood that the most important variables controlling snowmelt are radiative fluxes, energy advection, turbulent heat fluxes, and the temperature departures that synthesize these (Groisman et al. 1994; Zhang et al. 1997; Aizen et al. 2000; Ohmura 2001). Furthermore, the time period in which snow can begin melting is generally controlled by insolation, whereas interannual variation in this date can mostly be attributed to variability in downwelling longwave radiation, which is largely a function of cloud cover variations and heat and moisture transport changing the mean atmospheric thickness (Zhang et al. 2001). Regionally, radiative fluxes have been found to play a larger role in melt energy at high latitudes, with advective energy and resultant sensible heat fluxes contributing more to melt at lower latitudes and earlier in the season (Mioduszewski et al. 2014; Ohmura 2001; Leathers and Robinson 1997; Zhang et al. 1996, 1997). Turbulent fluxes can be important in the lower latitudes by enhancing or counteracting radiative fluxes (Male and Granger 1981) but have not been found to be significant drivers of snowmelt on a large scale (Zhang et al. 1997; Shi et al. 2013).
Synoptic conditions that generate the patterns that control energy advection have been studied on a regional and hemispheric scale (Aizen et al. 2000; Bamzai 2003; Vicente-Serrano et al. 2007; Ueda et al. 2003; Tedesco et al. 2009; Shi et al. 2011, 2013). There has been some success correlating winter snow conditions and subsequent melt season timing with teleconnections such as the Arctic Oscillation and Pacific–North American pattern, with up to 50% of the variance in these conditions explained in some regions (Tedesco et al. 2009; Bamzai 2003; Brown and Goodison 1996). Spring snow conditions also correlate with geopotential heights and modes of atmospheric circulation (Vicente-Serrano et al. 2007; Stone et al. 2002; Bao et al. 2011), but correlation strengths again depend on the region as well as methodological considerations such as time lags and spatial and temporal averaging. For example, using monthly averages not only removes some of the signal but also introduces mixed environmental conditions by incorporating both snow-covered (melting and frozen) and snow-free surfaces into the analysis, particularly one that covers many melt seasons.
A detailed analysis of energy balance terms is another method for understanding the snowmelt process. While synoptic conditions generally control these terms, atmospheric teleconnections are only a metric for certain modes of regional circulation, which are then manifested through the surface energy balance. Small- and point-scale studies of snowmelt attribution often have the advantage of utilizing energy balance data at a high temporal resolution in detailed analysis (e.g., Sicart et al. 2006; Stone et al. 2002; Pomeroy et al. 2006; Marsh et al. 2010) but can be difficult to generalize beyond the unique geography of the study location and are often limited to one melt season. Larger-scale energy balance studies have had some success attributing spring snowmelt while operating within a different set of constraints. Iijima et al. (2007) concluded that the advection of heat and moisture played the greatest role in eastern Siberian snow ablation, while Ueda et al. (2003) found that early snowmelt years on the eastern European plain are associated with anomalous low-level warm air advection.
Here, we use atmospheric reanalysis data to assess trends in the melt onset date (MOD) and its drivers, and to attribute MOD variability in distinctive regions. Regions of homogeneous snowmelt are identified using principal component analysis (PCA) and much of the analysis proceeds using these regions. Thirty-four-year trends in MOD and melt forcing variables are obtained using the Mann–Kendall trend test and Sen’s slope to assess regional MOD and its drivers. An attribution analysis is further performed using principal component regression, correlations between component loadings and mean surface atmospheric pressure, and a more detailed analysis of the influence of cloud cover and energy advection on MOD. Maintaining the focus on the surface energy balance allows for a unique and improved understanding of what type of forcing is necessary to initiate the snowmelt season, while combining these multiscaled approaches over a relatively long time period places results in the context of the observed hemispheric changes in snowmelt onset.
The study area includes all land between 55° and 75°N (hereafter simply Northern Hemisphere). The study area is divided into two regions for display purposes: the North Eurasian region is denoted E and is bounded by 40–180°E, while the North American region is denoted N and bounded by 60–170°W. Snow cover characteristics across this large area tend to be almost exclusively tundra and taiga as defined by Sturm et al. (1995). All of the Arctic coast, particularly Canada north and west of Hudson Bay, is characterized by bare or shrubby tundra, which transitions into generally taiga of varying canopy height and density in the remainder of this area (e.g., Bicheron et al. 2008). This area is further divided into subregions of homogeneous snowmelt onset for analysis, the methodology of which is described in section 2c (Fig. 1).
2. Data and methods
a. MERRA data
Near-surface energy balance and atmospheric variables, hereafter forcing variables, were obtained from the National Aeronautics and Space Administration (NASA) Modern-Era Retrospective Analysis for Research and Applications (MERRA) products (Bosilovich et al. 2011; Cullather and Bosilovich 2011, 2012; Rienecker et al. 2011). MERRA is run on a 1/2° latitude by 2/3° longitude global grid with 72 hybrid-sigma vertical levels to produce analyses at 6-h intervals covering the modern satellite era from 1979 to present. Data were obtained for a 122-day period from 1 March to 30 June, capturing the melt season from 1979 to 2012. The earliest mean MOD is in mid to late March, so there are a few grid cells during the earliest melt years when the MOD is earlier, but the effect of this on the analysis is insignificant. Different MERRA variables were used for different analyses depending on what was most appropriate, and include net and incoming shortwave (SW) and longwave (LW) radiation, LW and SW cloud radiative effect (CRE), sensible heat fluxes (positive fluxes moving from the snow surface to the atmosphere), vertically integrated atmospheric energy convergence (EC; a measure of energy advection), 2-m daily mean and maximum temperature, diurnal temperature range, 850-hPa temperature and specific humidity, and 1000–500-hPa thicknesses. These were aggregated from hourly to daily, although daily temperatures were derived from hourly 2-m temperature data. Statistical significance of trends at α = 0.05 in forcing variables as well as the date of melt onset were analyzed using a Mann–Kendall trend test (Mann 1945; Kendall 1975), with the associated Sen’s slope applied to the data to obtain a change over the 34-yr period (Helsel and Hirsch 2002).
MERRA has been evaluated extensively since its release (Cullather and Bosilovich 2012; Reichle et al. 2011; Robertson et al. 2011; Kennedy et al. 2011; Zib et al. 2012) but is subject to many of the limitations of other reanalysis systems. Because most of the data assimilated into MERRA are from remote sensing products, calibration issues across different platforms in the last 35 yr potentially introduce some bias into the outputs. These tend to be strongest over lower latitudes and over the ocean (Bosilovich et al. 2011; Robertson et al. 2011), with the greatest influence in the Arctic likely being negative biases in water vapor and downward LW radiation (Bosilovich et al. 2011; Kennedy et al. 2011). However, many of MERRA’s energy balance and advective terms used in this study have demonstrated some of the lowest biases among reanalysis products in comparison and validation studies (Cullather and Bosilovich 2011; Zib et al. 2012; Lindsay et al. 2014).
b. Melt onset data
The annual melt onset date data were obtained from Wang et al. (2013). Here, we briefly introduce the satellite data and the melt detection algorithm. Satellite-borne microwave sensors are effective tools for examining changes in snowmelt dynamics over the Arctic due to their high sensitivity to liquid water in snow and the general absence of cloud cover issues faced in visible imagery analysis (e.g., Glen and Paren 1975; Zwally and Gloersen 1977). Melt onset was determined using an algorithm based on temporal variations in the differences of the brightness temperature (Tb) between 19 and 37 GHz (vertical polarization) from passive microwave satellite data (Wang et al. 2013). A continuous time series of daily Tb was obtained from the Scanning Multichannel Microwave Radiometer (SMMR; 1979–87), SSM/I (1987–2008), and the Special Sensor Microwave Imager/Sounder (SSMIS; 2009–12) mapped to the 25 km EASE-Grid available at the National Snow and Ice Data Center in Boulder, Colorado (Armstrong et al. 2013; Knowles et al. 2002). Sensor cross calibration was performed by applying adjustment coefficients derived in previous studies (Abdalati et al. 1995; Jezek et al. 1993; Stroeve et al. 1998). Since our melt detection algorithm only uses the relative change in the temporal variations in Tb, slight offsets in Tb between sensors should not affect algorithm performance. Nevertheless, the differences in the SMMR and the SSMI/I and SSMIS sensors may still introduce some uncertainties in the long-term trends of MOD, and changes in sensor calibration could potentially have caused an apparent but not real change in MOD. As explained in Tedesco et al. (2009), this uncertainty cannot be quantified, but these data were calibrated as well as possible by Wang et al. (2013). The gaps in the SMMR data (due to a narrower swath and availability of every other day) and SSM/I and SSMIS data are filled by linear interpolation from adjacent days.
The melt detection algorithm is capable of distinguishing early periodic snowmelt onset from the main seasonal melt onset (Wang et al. 2008, 2013). Multiple melt events were determined based on temporal variations in the differences of the brightness temperature (Tb) between the vertical (V) polarized signals at 19 and 37 GHz (TbD = Tb19V – Tb37V) from the passive microwave satellite measurements. The MOD was detected if the difference in daily TbD and the previous 3-day average (M) was greater than a threshold (TH1 = 0.35 × M) for four or more consecutive days. An iterative approach was used for melt end date (MED) detection: melt end was detected as the first day when daily TbD was less than a threshold (TH2 = mean TbD in July + 7K) for at least 28 continuous days. If these conditions were not met at a given grid cell, the number of days was reduced from 28 to 21, and then to 14 if necessary. Melt duration was the number of days between MOD and MED for each melt event, with the main melt event identified as the event with the longest melt duration. More details on the determination of the melt detection thresholds can be found in Wang et al. (2013).
The melt detection algorithm was evaluated using observations at weather stations across the pan-Arctic region (Wang et al. 2008, 2013). The results showed that the primary MOD was associated with the early stage of the final ablation of the snowpack when the snowpack was wet but still fully covering the ground. The detected MODs corresponded to a clear shift in the statistical distribution of mean daily air temperatures from largely below freezing to above freezing. Since melt characteristics and melt season timing over permanent snow and ice are different from seasonal snow cover over land, a land ice mask was used to mask out those areas in our analysis (Brown et al. 2001). Because of the uncertainty of microwave measurements in complex alpine terrains (e.g., Tong et al. 2010), the performance of the melt detection algorithm in mountainous areas may have a larger uncertainty and needs to be tested further.
For melt attribution analysis, this MOD dataset is regridded to the MERRA model grid using a patch recovery interpolation method (Zienkiewicz and Zhu 1992). This interpolation method is a type of finite-element method that operates by determining (recovering) the derivatives of the finite elements at each node of the generated data mesh. Although computationally more expensive than more common piecewise polynomial interpolation methods, the “superconvergent” property of patch recovery allows for less resulting interpolation error because the solution error converges faster than expected based on the maximum error in the domain of the boundary value problem.
c. Principal component analysis of variability in the date of melt onset
Spatial clusters with similar MOD were identified with principal component analysis (PCA) in S mode, which is used to identify spatial clusters in time series (Richman 1986). Here, PCA in S mode was used on the standardized and mean-centered melt onset field to isolate regions of homogeneous variance in MOD. The principal components of PCA are eigenvectors, which are the linear combinations of the original variables but weighted by their contribution to explaining the variance in each dimension. These components were then linearly transformed, or rotated, using the varimax method (Kaiser 1958), implemented using MATLAB’s “rotatefactors” function. This is an important step for such an application of PCA because rotated components are then less domain-dependent and more localized in space, which helps to draw out physically meaningful clusters (Richman 1986).
To identify regions with similar MOD, each of the resulting 34 components was correlated with the melt onset anomaly field. In 16 of the 18 first components, one region of uniform correlation coefficients with R2 > 0.2 could be identified (Fig. 1). For this function of PCA, the strength of the correlation is not nearly as important as the spatial orientation of correlation coefficients. Components 19 to 34 explain less than 20% of the variance in the melt anomaly field as the components begin to degenerate into primarily statistical noise (i.e., are not physically meaningful; Hwang and Nettleton 2003). Not all Northern Hemisphere land is incorporated into the subregion given the nature of PCA, but we sacrifice some areal coverage for accuracy in the development of targeted regions. For example, regions with a strong gradient in MOD due to topography particularly near the coast will not be identified by the PCA algorithm, such as across much of the Rockies and in Greenland and Scandinavia. This limitation was part of the justification for not incorporating the longitudes of the latter in the broader analysis.
In the 16 regions identified, melt onset timing covaries among the region’s grid cells reasonably well through the 34-yr period. The MOD does not necessarily occur on the same date within each region every year, but it is likely that these regions respond to the same atmospheric forcing much of the time around the MOD. In many of these regions, there is considerable similarity in topography and land cover. For example, region 5N aligns well with the Yukon River Valley, 6N covers the tundra of northeastern Canada and part of the Canadian Archipelago, and 6E represents a significant subset of the Russian Plain. However, the Urals nearly cut region 4E in half, and region 8E overlays part of the Russian Far East with topography dominated by a complex system of mountain ranges, so topography is not always uniform. The extent to which land cover types are homogeneous throughout regions is more likely a function of climate rather than any influence land cover has on MOD, which would manifest in the PCA output. Additionally, dominant snowpack (Sturm et al. 1995) and land cover types are not nearly as important for snowmelt onset as they are for the melt dynamics of the ensuing melt season.
d. R-mode principal component analysis and multiple linear regression
Multiple linear regression analysis was performed with eight forcing variables encompassing a range of potential melt drivers: downwelling LW, incoming SW, sensible heat fluxes, LW CRE, 2-m temperature, 850-hPa specific humidity, 1000–500-hPa thicknesses, and EC. Anomalies were calculated relative to a fixed day of the year (DOY), that is, relative to the 34-yr mean of the variable for a given DOY. Then, the variable anomalies were averaged over a 1–4-day period prior to MOD such that the averaging period maximized the correlation between the 34-yr time series of time-averaged anomalies and melt onset anomalies. Anomalies of net LW and SW could not be used because of their strong dependence on surface conditions (i.e., snow cover), which strongly influences the relationship between anomalies and MOD timing, so incoming LW and SW were used instead.
Before regression was conducted, R-mode PCA was done on the eight variables. This type of PCA is the extraction of features from a set of variables such that the dataset is represented by a smaller number of new, orthogonal components (Richman 1986), rather than the S mode, which finds clusters in a spatial field of time series. Each component often comprises only one or a few highly correlated variables in the dataset, allowing physical interpretation while reducing dataset dimensionality and multicollinearity. They can then be used to replace the original variables in an analysis such as multiple linear regression, which is done here using MOD anomalies as the dependent variable. We seek to use these distinct “melt schemes” developed by the PCA to explain MOD anomalies through regression. For example, increased water vapor in the atmosphere does not alone drive melt onset in a particular year, but it can do so in coincidence with increased atmospheric thicknesses, surface temperature, LW radiation, and potentially low cloud cover.
Identifying the variables that comprise each PC is done by correlating each component’s loading pattern with the field of variables. The loading pattern is itself a correlation matrix, obtained by correlating the original data with each principal component, such that each component has a loading pattern. This second correlation results in a vector of correlation coefficients for each component, but determining the correlation cutoff value for which a variable should be considered to comprise that component is not trivial. This variable selection has been studied extensively (e.g., Cadima and Jolliffe 1995; Al-Kandari and Jolliffe 2001), and a universal cutoff value is not recommended (Jolliffe 2002). Therefore, we apply an empirical method to determine this cutoff loading by using biserial correlation analysis following Richman and Gong (1999), which locates the point at which the PC loading pattern agrees best with the corresponding signal in the original data. In this analysis, the biserial correlation peaks near a component loading value of 0.45, meaning that variables within each component with a correlation coefficient higher than 0.45 are considered to comprise that component.
A primary advantage of R-mode PCA is the reduction of variables made possible by constructing much of the variance in the data into the first few components. Several methods for component selection criteria in PCA exist [see, e.g., an overview in Jolliffe (2002)]; here, we selected the first components that explain more than 90% of the dataset variance. This resulted in retention of the first four components, and rejection of trailing components. In addition to containing redundancy and degenerating into noise, the trailing components have very small variances that would generate significant instability (i.e., large confidence intervals) in linear regression coefficients if they were to be used. Regression analysis proceeded on the retained components using a backward selection linear model, discarding resultant regression coefficients not significant at α = 0.05. Because both dependent and independent variables were standardized to run the PCA, the regression coefficients are in standardized units, and can be interpreted as the MOD response (standard deviation of days) to a one standard deviation change in the given component with the remaining components held constant.
e. Surface atmospheric pressure analysis
Synoptic conditions at MOD were identified by correlating the S-mode PCA loading pattern corresponding to each region with mean sea level pressure (SLP) averaged between the dates of earliest and latest melt in the given region. This reveals the spatial configuration of SLP during the range of time when snow begins to melt in each region. A positive correlation indicates the tendency for positive SLP anomalies during late melts and/or negative SLP anomalies during early melts, with the opposite true for negative correlations. Wind direction can then be inferred based on the orientation of SLP centers relative to each other. As such, this analysis also provides a physical mechanism for regional variations in EC that are explored further as well. However, this type of analysis does not account for climatologically favored, semipermanent systems such as the Siberian high and Aleutian low and their existing effect on wind direction.
a. MOD and forcing variable trends and variability
The mean MOD varies from March across the boreal forest region to early June in the Arctic tundra, with a strong inverse relationship with longitude (in addition to the first-order relationship with latitude) evident on both continents (Fig. 2a). MOD variability is largest across southern North America, northwestern Eurasia, and coastal locations. In those regions, the standard deviation exceeds two weeks (Fig. 2b), which is roughly consistent with the regions of earliest MOD (Fig. 2a). This pattern exhibits larger variation over locations that average earlier MODs, and about half as much variation over northern North America and eastern Eurasia.
There are marked differences in spatial patterns of the mean of forcing variables in the 3-day time period prior to and including MOD (Fig. 3). Most spatial patterns in these means are similar to the pattern in MOD, with the exception of LW CRE and, to a lesser extent, specific humidity and atmospheric thicknesses. More surface energy, primarily from radiative fluxes (Figs. 3a,b) and warmer temperatures (Fig. 3h), is available in the northernmost locations where snow begins to melt in May and June. This is primarily due to greater insolation closer to the solstice, with an average of over 150 W m−2 more net SW radiation compared to the southern part of the Eurasian and North American regions where melt starts earlier. Additionally, more moisture in the atmosphere on average in these northern regions, shown in 850-hPa specific humidity (Fig. 3e), combines with higher atmospheric thicknesses to generate levels of net LW radiation comparable to other regions, which would otherwise be lower due to less energy converging into the atmospheric column. This results in higher 2-m temperatures but a lower diurnal temperature range due to similar daily maximum temperatures near freezing.
A larger percentage of energy is derived from outside the region than from local radiative fluxes in southern and western regions of North America and Eurasia where snow begins melting earlier in March and April. This is shown in higher levels of EC farther south and west in these regions (Fig. 3g) but with lower atmospheric thicknesses (Fig. 3d) and specific humidity (Fig. 3e) leading to lower mean 2-m temperatures (Fig. 3h). While spatial differences in EC exceed 300 W m−2, only a fraction is available as part of the surface energy balance to melt snow. Some of this energy is radiated out to space or is temporarily stored in the air column, while the remainder heats the column and is thermally radiated to the surface. Finally, regional differences in LW radiation generated by clouds are under 40 W m−2 (Fig. 3f), although interannual differences can exceed 100 W m−2 given that years with relatively clear and overcast conditions are averaged together.
Most of the strongest trends toward earlier MOD are located across northern regions and are limited to regions 4E–10E, where melt begins 10–15 days sooner (Table 1). Across much of the remainder of Eurasia and northern North America, statistically significant (α = 0.05) trends suggest that MOD occurs 7–10 days earlier. The southern half of the North American region has no significant MOD trend, and no significant trend toward later MOD was found in any region. Significant trends in forcing variables exist mostly in regions where there are strong trends in MOD (Table 1). Because of the strong seasonal cycle of forcing variables, earlier MOD occurs when most forcings (particularly net LW and SW, specific humidity, 1000–500-hPa thickness, and 2-m temperature) have smaller values than at the date of mean melt onset (Fig. 4). The only exception is EC, which is higher earlier in the year when the atmosphere acts as more of a sink for energy. Thus an earlier shift in MOD results in an increasing trend for EC. However, it is unclear whether increased EC is driving earlier snowmelt, or whether it only appears larger because it tends to be earlier in the spring; this is investigated further below. Additionally, the lack of a trend in LW radiation in all regions suggests that these opposing trends between EC and other correlated terms such as atmospheric thicknesses and specific humidity that influence LW are negating each other.
b. Principal component regression
Results from the R-mode PCA show that the first principal component comprises a combination of LW radiation, specific humidity, 1000–500-hPa thicknesses, and 2-m temperature in most regions, explaining 45%–55% of the variance in the data field (Table 2). Although more variance is contained in the first component, this only indicates that these variables comprise the most spatially coherent covariance in the data with the trailing components being constrained orthogonally, and they are not necessarily the most important predictor of MOD anomalies. The second component typically describes cloud cover, with LW CRE, incoming SW (of the opposite sign), and sometimes specific humidity exceeding the 0.45 correlation coefficient truncation level. LW radiation often appears significant in both components because it is highly correlated with water vapor, atmospheric thicknesses, and 2-m temperature, but also has a strong positive and negative correlation with LW CRE and SW radiation, respectively. The third and fourth components are both largely a combination of only sensible fluxes and EC.
The results of the regression show that all regression coefficients for PC 1 are statistically significant, and are consistently some of the largest, with absolute values ranging from 0.33 to 0.47 (Table 3). Regression coefficients for PC 2 are generally smaller, but with a positive sign, indicating that anomalously early (late) MODs are correlated with negative (positive) anomalies in cloud cover in most regions. The composition of the first and second components in regions 9E, 10E, 6N, and 3N is mixed between these two melt schemes (Table 3, bold); that is, both components exhibit a significant correlation with both sets of forcing variables that are otherwise split between these first two components (Table 2). Therefore, it is not clear how to interpret those regression coefficients, and this ambiguity also lowers some of the coefficients in PC 1 and 2. Only one component of PC 3 and 4 is typically significant in the regression, so the regression coefficient that is significant is shown regardless of whether it represents PC 3 or 4. Regression results for this component indicate that positive (negative) EC anomalies explain some early (late) melt variability primarily in regions with earlier mean melt dates, although northern North America is included as well.
c. Analysis of energy convergence
The influence of EC as a melt driver was explored further by region on account of the importance of advection in providing melt energy. EC time series in the week before and after MOD have a very distinct pattern where EC nearly always rises sharply to a peak just before MOD before returning to premelt levels (Fig. 5a). The maximum increase to its peak in the 3- or 4- day time period before MOD was at least 300 W m−2 on average in all regions, with the exception of regions 2E, 4E, and 1N where it increased at between 500–600 W m−2 (Fig. 5b). Similarly, EC peaks at over 300 W m−2 in all regions and achieves a peak of over 400 W m−2 across some of the earliest melt regions. Regional differences in EC anomalies (relative to fixed DOY as described previously) show the same pattern (Fig. 5b), with larger anomalies in regions farther to the south in North America and south and southwest in Eurasia.
The time series of EC prior to MOD shows considerable variability in terms of standard deviation (Fig. 5b), but differing levels of variability regionally. In earlier melt regions, greater variability in the EC time series exists and manifests as the observed larger increase just prior to MOD (Fig. 5a) and results in larger EC anomalies. Farther north, EC is greater prior to MOD with a smaller increase, and it follows that EC anomalies are lower here. These results show that regions with an earlier MOD exhibit more variability in EC prior to MOD with a larger increase to a higher peak in the 3–4 days prior to MOD. As a result, EC anomalies are higher here, agreeing with regression results showing that EC anomalies explain the most variability in MOD in many of the same regions.
d. Surface pressure analysis
The correlation between SLP around MOD and S-mode loading patterns provides insight into the regional preference for certain synoptic conditions if they are very different between early and late melts. The spatial configuration of SLP favors flow from the south or southwest during years with early melt onset in many regions such as regions 2E–7E where EC has been shown to be a more significant melt driver (Fig. 6). This is demonstrated by positive correlations to the west or northwest of the given region and/or negative correlations to the southeast or south, indicating the preference for low and high pressure, respectively, during earlier melt years. The opposite interpretation could be constructed for late melt years, though early melt forcing tends to be more anomalous. Additionally, the climatological mean flow (not shown) is from the south and southwest during March and April across southern and western Eurasia, in the same regions as mentioned above, further supporting the greater influence of EC.
The synoptic pattern is unclear across the northernmost regions where there is either little correlation between SLP and MOD, as in regions 9E and 6N, or a more zonal flow as in regions 10E and 5N. Regions 2N and 3N in central Canada exhibit a very similar correlation pattern to each other. However, with a large area of anomalous lower SLP to the northwest during early melts (or higher SLP during late melts), the effect of this SLP configuration on the circulation is unclear; this is too far displaced to be a mode of the Pacific–North America pattern. The general circulation does not appear to favor more energy advection in early melt years (or less in late melt years) in these regions. This supports the conclusion that the temperature increase to the freezing point at MOD occurs through means other than EC.
e. Analysis of cloud cover
Because of its ability to significantly alter the surface radiation balance, the influence of clouds was examined in more depth by evaluating the change in LW CRE prior to MOD. LW CRE was found to decrease by at least 10 W m−2 approximately twice as often as it increases past that threshold before MOD, with a corresponding increase in net SW radiation (Fig. 7). The exceptions to this are regions 1E, 6N, and 4N where LW CRE increases about as often as it decreases, although these three regions are generally dissimilar otherwise. While the magnitude of change is not directly calculated, this gives an interannual cross section of the typical cloud cover conditions rather than averaging LW CRE over the 34-yr period. Furthermore, given that the increase in net SW radiation occurs slightly more often than a decrease in LW CRE, it is reasonable to assume that the magnitude of the corresponding increase in SW radiation offsets or exceeds the loss in LW from clouds more often than not.
Multiple linear regression and PCA analysis show that the strongest predictor of early (late) MOD anomalies are positive (negative) anomalies of atmospheric thicknesses, water vapor, and consequently LW radiation and 2-m temperature, and that regional differences in this regression coefficient are small. These four forcing variables comprise PC 1 in all regions except 9–10E, 3N, and 6N, and as such have the strongest covariance between 1979 and 2012. EC has less coherent anomalies and higher-frequency variability, and exhibits insignificant correlation with the components until PC3 and PC4 when 65%–75% of the data’s variance has been accounted for. Positive EC anomalies and sensible heat flux are associated with early MOD anomalies, particularly in southern and western Eurasia. However, EC anomalies do not always correlate well with MOD anomalies, which makes the regression less useful for an analysis of its influence. Nonetheless, PC3–PC4 regression coefficients are largest through southern and western Eurasia as well as western North America where EC anomalies peak immediately prior to the MOD, whereas coefficients are not significant across northern Eurasia and much of North America.
While much of the analysis of regression coefficients is focused on explaining melt in early melt years, melt onset in late melt years is still initiated by some forcing. We hypothesize that this forcing is most likely to be the increased amount of insolation that more readily brings the temperature to the freezing point. Insolation represents a forcing that exists to initiate melt if EC does not do so first because it does not provide a mechanism by itself; its daily increase is minimal. But SW radiation adds over 150 W m−2 on average (more in those regions when melt begins later) after accounting for surface reflection, which is a significant contribution to the surface energy balance and the energy available for snowmelt. Although there is some increase in LW radiation throughout the season from higher geopotential heights and water vapor as the general circulation transitions into a warm season regime, the consistent increase in insolation tends to exceed this as shown in Fig. 3 (evidence is indirectly based on mean MOD) and Fig. 4.
Our conclusion that insolation is more important both in late melt years and in regions where melt begins later is consistent with other research, particularly in contrast with turbulent fluxes. It is generally agreed that turbulent fluxes have a minimal influence in the Arctic and sub-Arctic, with a greater influence in the midlatitudes and marine locations, while insolation has a much greater effect farther north and with later melts (Leathers and Robinson 1997; Zhang et al. 1997; Ohmura 2001; Liston and Hiemstra 2011; Shi et al. 2013). While we are unable to separate sensible fluxes from EC in the regression coefficients, the magnitude of average sensible fluxes is low, and of positive sign (directed away from the surface) across much of the study area (Fig. 3b). There is larger diurnal variation in this term, so a positive mean value for sensible flux is still likely directed toward the surface for part of the day, but its overall contribution to melt energy is small.
Cloud cover in all regions except 1E, 6N, and 4N decreases prior to MOD in at least two-thirds of the years, and this decrease is independent of the magnitude of the MOD onset anomaly. In this study we find that on average, the increase in insolation when clouds disappear compensates for or exceeds the loss of LW radiation, thus neutralizing the energy balance effect. Regression results showing more cloud cover correlating with late melts (Table 2) is more due to this decrease rather than a causal relationship with MOD anomalies, since LW CRE anomalies in late melt years tend to be higher relative to the same DOY in other years. The ability of clouds to substantially add LW energy to the surface and initiate melt is well documented (e.g., Zhang et al. 1997; Stone 1997; Stone et al. 2002) and is a factor in some years. However, this study shows that this likely occurs in a minority of years in most of our study regions, and the surface from an energy balance standpoint does not depend on a surge of LW radiation to begin the melt season.
We show that there are large regional differences in EC such that much of the southern part of the entire study area relies more on EC as a melt driver. In these locations, EC is greater, more variable, and more anomalous both in regions where melt begins earlier and in anomalously early melt years. Much of the magnitude in these anomalies is generated only a few days prior to MOD, as shown in Fig. 5. Furthermore, there is a trade-off between EC and radiative fluxes (primarily SW, as LW is also generated from moist static energy convergence into the column) in terms of what is present at MOD; radiative fluxes play less of role in melt in locations where EC is the dominant driver. There is also considerable zonal asymmetry in meridional energy transport, with poleward transport climatologically occurring at preferred longitudes such as western Siberia and Alaska, with even a net equatorward transport over northeastern Canada. This was calculated across the 70th parallel by Serreze et al. (2007) and attributed to differences in the mean longwave circulation patterns in the middle troposphere. This explains some of the inconsistency in Fig. 6 with other analyses; anomalies in SLP do not need to be strong to facilitate extensive energy transport across southern and western Eurasia where Figs. 3 and 5 show EC is strongest.
Our results show that regions that rely more on advective energy exhibit greater interannual variability in MOD due to the weekly to monthly duration of prevailing synoptic patterns, whereas regions dominated more by radiative energy should have a more predictable onset date. This is seen in Fig. 2b, with MOD standard deviations in western Eurasia and southern and eastern North America about twice those in northeastern North America and eastern Eurasia. Lack of advective energy in the latter regions is due to the gradual climatological shift in the circulation pattern through the spring as Rossby wavelengths lengthen and become less amplified, limiting opportunities for transient eddies to transport heat and moisture poleward.
Trends in MOD and many of the forcing variables differ among groups of regions. There are relatively strong trends toward earlier melt in regions 1E–7E with corresponding significant increases in EC and decreases in atmospheric thicknesses, specific humidity, and weaker trends downward in net SW radiation. Two-meter temperatures show trends downward by 2°–3°C, although some decrease is accounted for in a 1°–2°C increase in diurnal temperature range; 850-hPa temperatures also show a 2°–3°C trend downward, which precludes the possibility that warm air aloft is still advected into the region to generate LW radiation despite the colder surface temperatures. Therefore, we theorize that the increase in EC compensates for the lower remaining energy balance terms in these regions. This may be a primary forcing behind the earlier melt, but the effect would be the same: with less radiative energy available earlier in the season, the atmosphere naturally acts as more of a sink for energy, which manifests as a larger amount of EC. Regardless of the extent to which it is a cause or an effect, we have shown that there is more energy available to converge into the atmospheric column. Additionally, the strong signal for EC as a melt forcing across much of the southern and western Eurasian region is consistent with other studies showing strong correlations with southerly and southwesterly flow bringing warmer and more moist air into primarily western and central Siberia (Ueda et al. 2003; Vicente-Serrano et al. 2007)
In contrast to regions 1E–7E, forcings driving the strong trends toward earlier melt in regions 8E to 10E are not possible to identify. In these areas, insolation has a stronger control on MOD, and they are more insulated from anomalous energy transport from the south. EC is constant as MOD shifts to earlier dates, and other variables exhibit weak trends when melt begins in May and early June. Nonetheless, decreasing amounts of net SW radiation here and in many other regions where MOD is arriving earlier suggest there are some broad limits on how anomalously soon MOD could come in a given location. Even earlier spring melts would occur with climatologically less energy in the atmosphere and especially less from insolation, particularly close to the equinox when insolation changes rapidly. This would require an even larger fraction of energy to come from EC, and there is likely a theoretical limit to how much energy can be transported poleward, so this may limit how much earlier MOD can occur in some places.
In North America outside of region 6N, MOD trends are not significant, and there are few significant trends in forcing variables. Wang et al. (2013, their Fig. 2a) shows linear regression-based trends in the same MOD dataset over 33 years on a grid-cell basis over the Northern Hemisphere north of 50°N, which shows significant MOD trends in a sizable minority of Canada and Alaska. This demonstrates that MOD trend differences between continents are not as large as they may appear here due to region selection and regional pixel averaging. It has also been established that the trend in snow disappearance is greater than that in earlier MOD, so there may be a shortening of the melt season in some places, as suggested in Tedesco et al. (2009), who found a shortening of 0.57 days yr−1 across the Arctic due to the difference in trends between these two dates. The explanation for why there are few significant trends in forcing variables over North America is less clear, but changes in spring atmospheric circulation across northern Alaska (Stone et al. 2002) and northern Canada (Brown and Braaten 1998) have been documented. How these changes translate into warmer spring temperatures and earlier melt onset regionally is not a simple problem, and given the vastly different geography of these two continents there may be different mechanisms at work that require further research.
We have found that the predominant snowmelt drivers in the Northern Hemisphere can be broadly split into two categories: those processes involved in energy advection, and those in radiative fluxes. This is unsurprising given the range of latitude, land cover types, continentality, and consequently the more than 3-month difference in mean MOD. Further, melt drivers shift from advective toward radiative as the MOD occurs closer to the solstice with stronger insolation and weaker synoptic influences. This may account for the weaker variability in MOD in the northernmost regions where this MOD occurs later in May and June.
As snowmelt occurs earlier in the spring, the dominant melt drivers may shift more toward advective mechanisms depending on corresponding regional changes in atmospheric circulation. While insolation controls the time period when melt can begin, increased equatorward energy available for transport north will continue to force earlier melting across much of the Northern Hemisphere. Given the regional variation in MOD standard deviation, this shift would likely result in greater variability in each season’s MOD. Furthermore, reduction in snow cover at lower latitudes would further enhance energy advection northward, which may already be occurring here demonstrated by a shift toward advective drivers farther north. The impacts of this type of shift in spring snow cover on the rest of the cryosphere and boreal climate system remain uncertain.
An assessment of the energy and radiation balance variables driving melt onset across much of the Northern Hemisphere over a 34-yr period around the beginning of the snowmelt season has led to the following conclusions:
The melt onset date has advanced by 1–2 weeks in the 1979–2012 time period across much of the Northern Hemisphere, with the strongest trends in northern and western Eurasia. In these regions, an increasing amount of advective energy is replacing decreasing levels of many other forcing variables associated with an earlier MOD.
Large differences in forcing variables at MOD generally follow the spatial pattern of melt timing, with weaker radiative fluxes and colder mean 2-m temperatures where MOD occurs earlier, but a greater magnitude of energy convergence. This pattern is also found in trends: as MOD comes sooner, there is less energy in the system and more advected from elsewhere.
Evidence supports the hypothesis that the date of melt onset is less variable farther north and therefore relies less on EC and more on radiative fluxes (insolation as well as LW radiation generated from higher atmospheric thicknesses and water vapor). Conversely, interannual MOD variability is larger farther south in regions that receive more outside energy. This is consistent with MOD synoptic-scale pattern changes that can generate flow from the south and southwest and increased energy advection.
Linear regression and PCA analysis show that the strongest predictors of early (late) MOD anomalies are positive (negative) anomalies of atmospheric thicknesses, water vapor, and consequently downward LW radiation and 2-m temperature. Secondary to this, positive (negative) anomalies of EC are associated with early (late) MOD anomalies in many regions and regression coefficients are of similar magnitude to those of the first principal component across southern Eurasia in particular.
Regional differences in EC are substantial. Compared to other regions, southern and western Eurasia and southern North America stand out. Here, EC increases more to a higher peak prior to MOD, is more of an anomaly relative to DOY, and is more variable from day to day. This suggests that EC is a much more important melt driver in these regions than in northern Eurasia and North America.
There is little difference regionally in cloud cover at MOD with no significant change over the study period, and often a decrease in cloud cover just prior to MOD.
Snowmelt onset across the high latitudes of the terrestrial Northern Hemisphere appears to be controlled by relatively different processes. As snowmelt occurs earlier in the spring, the dominant melt drivers may shift more toward advective mechanisms depending on corresponding regional changes in atmospheric circulation. The uncertainty that remains in the observed and future variability of snowmelt requires further work to help resolve this. Given the limitations of reanalysis data and the difficulty in separating cause from effect in some of the processes involved in snowmelt, it is crucial to continue to address this important component of the changing cryosphere with a variety of methods and perspectives.
J. R. Mioduszewski was supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program, Grant NNX13AO44H. D. A. Robinson and J. R. Mioduszewski were also supported by NASA MEaSUREs Award NNX08AP34A. We thank three anonymous reviewers for their constructive suggestions.