Recent ice cap snowmelt in Russian High Arctic and anti-correlation with late summer sea ice extent

Glacier surface melt dynamics throughout Novaya Zemlya (NovZ) and Severnaya Zemlya (SevZ) serve as a good indicator of ice mass ablation and regional climate change in the Russian High Arctic. Here we report trends of surface melt onset date (MOD) and total melt days (TMD) by combining multiple resolution-enhanced active and passive microwave satellite datasets and analyze the TMD correlations with local temperature and regional sea ice extent. The glacier surface snowpack on SevZ melted significantly earlier (−7.3 days/decade) from 1992 to 2012 and significantly longer (7.7 days/decade) from 1995 to 2011. NovZ experienced large interannual variability in MOD, but its annual mean TMD increased. The snowpack melt on NovZ is more sensitive to temperature fluctuations than SevZ in recent decades. After ruling out the regional temperature influence using partial correlation analysis, the TMD on both archipelagoes is statistically anti-correlated with regional late summer sea ice extent, linking land ice snowmelt dynamics to regional sea ice extent variations.


Introduction
Mass wastage of glaciers and ice caps significantly contributes to the present-day global eustatic sea level rise, and is strongly influenced by regional climatic changes (Bahr and Radić 2012, Jacob et al 2012, Kaser et al 2006, Meier et al 2007. Air temperature rise in the Russian High Arctic is outpacing many other Arctic regions with an average magnitude of 1°C to 3°C since the mid-1950s (IASC 2010, IPCC 2007, Kalnay et al 1996, exacerbating the vulnerability of the small and low-elevation glaciers on Novaya Zemlya and Severnaya Zemlya to mass loss (figure 1). Surface melting plays an important role in the interannual and decadal variability of surface mass balance of Russian High Arctic glaciers (Dowdeswell and Hagen 2004, Bassford et al 2006, Sharp and Wang 2009, Zeeberg and Forman 2001. Snowmelt also significantly impacts regional climate by lowering surface albedo, increasing the absorption of solar radiation and amplifying the positive snow-temperature feedback (Karner et al 2013, Mätzler 1994, Sharp and Wang 2009). However, due to the harsh environment, field studies on broad-scale surface-melt dynamics on Novaya Zemlya and Severnaya Zemlya are spatially or temporally limited. As a powerful alternative, remote sensing technology has been applied in detecting melt processes in this region; however singlesatellite measurements usually cover fewer than ten years because of short-term operation Wang 2009, Smith et al 2003), which constrains the understanding of long-term surface-melt response to recent rapid warming in Russian High Arctic. By combining multiple resolutionenhanced remote sensing datasets in the present research, we are able to present a continuous decadal measurement  and trend analysis of snowmelt onset date (MOD) and total melt days (TMD) for the ice caps on Novaya Zemlya and Severnaya Zemlya.
Recent sea ice shrinkage has been identified as a potential forcing on glacier surface mass balance through its influence on precipitation patterns (e.g. Carr et al 2013). Bamber et al (2004) suggested that the mass growth anomaly on the central Austfonna ice cap on Svalbard during 1996-2002 was associated with a simultaneous precipitation increase, which was most likely caused by an enhanced moisture transport resulting from declining perennial sea ice cover in the Barents Sea. Model experiments also suggested that decreased sea ice concentrations in the Southern Ocean would cause higher surface temperatures over greater open water areas and enhanced evaporation from the surface, which can contribute to greater precipitation along the coast of Antarctic (Weatherly 2004). A recent study reported positive trends in precipitable water over much of the Arctic, possibly as a combined result of negative anomalies in late summer sea ice extent and positive anomalies in air and sea surface temperatures (Serreze et al 2012). Furthermore, Arctic sea ice was recognized as a contributing driver of recent northern European summer precipitation increase by altering the largescale atmospheric circulation (Francis 2013, Screen 2013. Sea ice also influences glacier and ice cap surface melting. Koerner (2005) proposed that increasing open water in Jones Sound might form melt-suppressing fog and result in the slower ablation-rate increase between 1961 and 2003 at lower altitudes of Sverdrup Glacier, Canada. However, Rennermalm et al (2009) suggested that reduced offshore sea ice concentration can warm the ocean mixed layer and promote onshore advection of latent heat and sensible heat, further increasing snowmelt extent on Greenland ice sheet. Regions shown from left to right in ordered stipple, simple hatch and crosshatch patterns are Barents Sea, Kara Sea and Laptev Sea, whose sea ice extents are correlated with snowpack melt on Severnaya Zemlya and Novaya Zemlya. Sea ice sectors come from the Arctic ROOS web-portal www.arctic-roos.org. An automatic weather station (AWS) was maintained during an ice core drilling campaign from 1999 to 2000 on Akademii Nauk ice cap in Severnaya Zemlya at 80°31'N, 94°49'E (Opel et al 2009 The wide-spread Greenland ice sheet melt in July 2012 followed the major acceleration of Arctic sea ice retreat in June (Hanna et al 2013), and Screen (2013) suggested that reduced sea ice in summer favors enhanced ridging formation over western N. Atlantic, which might have contributed to the expanded 2012 Greenland surface melt event (Francis 2013, Nghiem et al 2012. Outside Greenland, Rotschky et al (2011) indicated the possibility of sea ice conditions as a driving force for glacier melt in Svalbard and found that, for some abnormal years, early snowmelt and long melt duration were related to local, anomalously low sea ice extent. Since 1978, sea ice reduction in the Russian High Arctic is among the fastest-declining regions in the Arctic (Cavalieri and Parkinson 2012) and Sharp and Wang (2009) pointed out that the heat released from an increase in surrounding open ocean waters might prolong the melt season of Novaya Zemlya and Severnaya Zemlya, but they lack long-term observations to establish the association conclusively. All of these observations motivate us to test the validity of the seemly apparent but evidence-lacking correlation between sea ice extent and glacier melting in this region. Put simply, we investigate the possibility of sea ice extent as an independent factor contributing to the melting of ice caps, which is particularly important for predicting how Russian High Arctic glacier surface mass balance will change if surrounding sea ice extent continues to decrease in the future. Since sea ice extent and glacier melting can co-respond to the regional temperature increase, partial correlation analysis is employed here to remove the large-scale climate influence on both variables. We explore the possible correlation using our retrieved total melt days and publicly available regional sea ice extent and reanalysis temperature data.

Study area
The present study focuses on the heavily glaciated archipelagoes of Novaya Zemlya (NovZ) and Severnaya Zemlya (SevZ) in the Russian High Arctic (RHA) (70°∼ 83°N, 45°∼ 110°E) (figure 1). Franz Josef Land is excluded from the present research because passive microwave signals of snowmelt over the islands contain significant noise from ocean water. SevZ and NovZ span a wide range of longitude and encompass more than 40 000 km 2 of glaciated area combined (IASC 2010, Moholdt et al 2012. Polythermal glaciers are dominant on the two archipelagoes (Glasser 2011). The annual mean air temperature ranges from −5°C to −15°C across this region with a July mean temperature slightly above 0°C in SevZ (Bassford et al 2006) and the northern tip of NovZ (Kotlyakov et al 2010). The climate is relatively mild and humid in southwest NovZ and gradually becomes drier and colder northeast into SevZ. NovZ is in the path of the warm North Atlantic Drift, and the prevailing westerlies promote warm air advection, bringing moist air masses from the Norwegian Sea and resulting in higher precipitation and an overall July mean temperature of 6.5°C on NovZ (Moholdt et al 2012, Zeeberg andForman 2001). Geographically, SevZ divides the Laptev Sea from the Kara Sea and NovZ is located in the middle of the Kara and Barents Seas ( figure 1). The surrounding open ocean waters serve as major moisture source for both archipelagoes (Opel et al 2009, Zeeberg and Forman 2001.

Data and methods
3.1. Data 3.1.1. Glacier area. We use the Russian Arctic outline in the Randolph Glacier Inventory (version 2.0), provided by Global Ice Measurements from Space (GLIMS) (Raup et al 2007), to delineate the glacier cover in the study area (figure 1). This version of the dataset was updated in June 2012 and has been frequently applied in recent glacier-related studies. To minimize mixed pixel effects from land and ocean (Bellerby et al 1998), we focused our analysis on ice cap interiors, covering approximately two thirds of the reported glaciated areas on each archipelago (see supplementary notes and figure S1 in the supporting information for more details on the narrowed-down study area).
3.1.2. In-situ and reanalysis temperature data. There are two meteorological stations with long-term records in the study regions, namely Golomjannyj Island (western SevZ: 79°3 6'N, 90°36'E) and Malye Karmakuly (southwestern NovZ: 72°24'N, 52°42'E) (figure 1), but to our knowledge, few temperature data were directly measured on ice in SevZ and NovZ during our study period. However, an automatic weather station (AWS) was run during an ice core drilling campaign in SevZ at 80°31'N, 94°49'E (700 m a.s.l., figure 1) and collected hourly data between May 1999 and May 2000 with data gaps in February and April 2000 due to power breakdowns (Opel et al 2009). The AWS recorded continuous daily temperature data covering the 1999 summer melt season and are employed here to relate to the microwave snowmelt signatures. Theoretically, the melting temperature of snowpack is determined from glacier surface energy balance; however, these AWS data do not have all necessary measurements to solve the energy balance equations. Therefore, the present work uses the above ground 2.5 meter air temperature data and assumes that melting occurs when daily maximum temperature exceeds 0°C following earlier similar snowmelt studies (e.g. Smith et al 2003, Trusel et al 2012. In order to put the surface melt dynamics into a climatic perspective, we determined large-scale temperature variations on NovZ and SevZ from NCEP-NCAR Reanalysis-1 data (Kalnay et al 1996). This dataset is used following earlier published ice surface snowmelt studies, such as ice caps in the Eurasian Arctic including NovZ as well as SevZ (Rotschky et al 2011, Sharp andWang 2009), Canadian Arctic (Wang et al 2005) and Greenland ice sheet (Wang et al 2007). We averaged the temperatures of reanalysis grid cells that cover the glacier areas and use the 'free air' 850 hPa geopotential height instead of the 2 m height in order to limit the sea surface temperature influence over the archipelagoes Environ. Res. Lett. 9 (2014) 045009 M Zhao et al (Moholdt et al 2012). In addition, temperature at 850 hPa geopotential height is found to be more correlated to glacier melt dynamics than temperatures at other levels in the Russian High Arctic (Sharp and Wang 2009). We also made a comparison between NCEP-NCAR and ERA-Interim (Dee et al 2011) reanalysis temperatures at 850 hPa geopotential height to discuss the confidence level of NCEP-NCAR reanalysis temperatures in our study region. The grid size is about 2.5°× 2.5°for NCEP-NCAR and 0.75°× 0.75°for ERA-Interim.
3.1.3. Resolution-enhanced passive and active microwave datasets. Visible band satellite images are limited in the RHA because of clouds and darkness (Marshall et al 1994); therefore, microwave data, which are highly sensitive to snow wetness, lack cloud contamination and are independent of solar illumination (Ashcraft and Long 2006), are ideal alternatives. There are two major types of microwave remote sensing sensors, namely active and passive. Active sensors transmit microwave energy to earth surface and measure the reflected energy as backscatter coefficient (σ 0 , unit: dB), which is a function of surface roughness, wetness, land type, polarization, frequency, topographic factors and sensor configurations (Hinse et al 1988). Instead of sending energy, passive microwave sensors directly collect earth surface self-emission as brightness temperature (T b , unit: K), which is approximated as a product of surface emissivity and physical temperature. Passive microwave sensors in the present study consist of Special Sensor Microwave Imager (SSM/I) and Advanced Microwave Scanning Radiometer for EOS (AMSR-E). K aband (∼37 GHz) vertically polarized T b channels are selected due to better correlations with land surface temperature and sensitivity to snow moisture (Apgar et al 2007, Kim et al 2011, Ramage and Isacks 2002. Active microwave sensors include the European Remote Sensing (ERS) Advanced Microwave Instrument (AMI), SeaWinds on QuikSCAT (QSCAT) and Advanced Scatterometer (ASCAT) onboard satellite Metop-A.
The original resolutions of passive and active microwave data are tens of kilometers and are challenging to use for surface melt studies on small ice caps such as in the NovZ and SevZ archipelagoes. The Microwave Earth Remote Sensing Laboratory at Brigham Young University (MERSL-BYU) resampled the original passive and active microwave data to enhanced (higher) resolutions using the Scatterometer Image Reconstruction (SIR) technique, which has been discussed and validated in Long and Daum (1998) as well as in Early and Long (2001). SIR technique is based on a multivariate form of block multiplicative algebraic reconstruction, which generates enhanced resolution gridded images from irregularly spaced measurements given sufficiently dense sampling (Long and Daum 1998). Resolution-enhanced passive microwave data have been successfully applied in detecting surface melting on Greenland and Antarctic ice sheets as well as monitoring sea ice dynamics (e.g. Agnew et al 2008, Ashcraft and Long 2005, Howell et al 2010, Kunz and Long 2006).
Resolution-enhanced active microwave data also have been widely used in global snowmelt studies (e.g. Rotschky et al 2011, Sharp and Wang 2009, Trusel et al 2012, Wang et al 2005, 2007. Therefore, the MERSL-BYU microwave datasets are utilized here to optimize for spatial resolution. The SSM/I F13 data from 1995 to 2007 and AMSR-E from 2003 to 2011 are utilized in the present research, and both of them are enhanced to 8.9 × 8.9 km 2 pixel grid resolution by MERSL-BYU. Active sensors' images from MERSL-BYU are used for a range of periods: ERS C-band (5.6 GHz) vertically polarized surface backscatter (σ vv 0 ) in 8.9 × 8.9 km 2 pixel grid from 1992 to 2000, QSCAT K u -band September sea ice extent is strongly correlated with mean summertime atmosphere temperature, and can be treated as a cumulative result of climatic forcing during the entire melt season (Stroeve et al 2011, Ogi andWallace 2007). Furthermore, monthly sea ice extent is auto-correlated and mean sea ice extent corresponding to TMD period is more sensitive to how much new sea ice was formed in the previous winter season, potentially introducing biases from earlier accumulation season. Therefore, the September sea ice extent is used here to assess the correlation between sea ice extent and total snowmelt days. The Laptev Sea, Kara Sea and Barents Sea sectorsʼ monthly mean September sea ice extents from 1995 to 2011 are provided separately by Nansen Environmental and Remote Sensing Center, Bergen, Norway through Arctic ROOS web-portal www.arctic-roos.org. The sea ice data are obtained from Scanning Multichannel Microwave Radiometer (SMMR) (1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987), SSM/I (1987SSM/I ( -2006, and Special Sensor Microwave Imager/Sounder (SSMIS) (2006present) using the Norwegian Remote Sensing Experiment (NORSEX) sea ice algorithm, which are based on radiative modeling of ice surface air temperatures and passive microwave channel ratios. The detailed algorithm description and validation are available in Svendsen et al (1983).

Methods
3.2.1. Passive microwave melt onset and total melt day algorithms. Wet snow has much higher microwave emissivity than dry snow, which results in an abrupt and distinct increase in a T b time series when surface snow changes from a dry state to a wet state (Stiles and Ulaby 1980). The minimum count of the bimodal T b distribution for perennially snow and ice covered surfaces is the threshold to separate the wet snow (higher T b ) from the dry snow (lower T b ) (Ramage and Isacks 2002). SSM/I 37 GHz vertically polarized brightness temperature (T b37V ) distribution of both islands from 1995 to 2007 indicates an annually consistent  figure 2(b)).
Melt onset for an individual pixel is usually the first date that the corresponding sensor's threshold (233 K for SSM/I and 252 K for AMSR-E) is exceeded for 3 of 5 consecutive days (e.g. Monahan and Ramage 2010). However, NovZ and SevZ glaciers tend to melt intermittently with melting and refreezing cycles throughout the entire ablation season, and abnormally short early melt events (Semmens et al 2013), usually occurring at the end of the accumulation season but preceding the ablation season, could affect the calculation of average snow melt onset date (MOD); therefore, the present work defines each archipelago's MOD as the first day when more than 90% of selected pixels melted (see supplementary notes and figure S5). Total melt days (TMD) are the sum of days that the corresponding sensor's threshold is exceeded, including short-lived events before or after the main melt season.
The difference between TMD for SSM/I and AMSR-E is calculated for the period of sensor overlap (2003)(2004)(2005)(2006)(2007) and shows a slight offset near-normal frequency distribution, indicating that SSM/I systematically detects about 7 more melt days than AMSR-E ( figure S4(a)). The different sensor sensitivity to snow wetness, melt thresholds, effective satellite overpass time and raw data footprint might account for the systematic error. Because of less variability in SSM/I melt threshold (see figure S2) and good agreement with in-situ observations ( figure 2(a), supplementary notes), we adjusted AMSR-E TMD values to SSM/I standard to account for the systematic difference (from 2008-2011) by adding 6.87 and 7.39 days to AMSR-E-derived TMD on SevZ and NovZ, respectively ( figure S4(a)). These combined records constitute a 17-year TMD record from 1995 to 2011 with a nominal resolution of 8.9 × 8.9 km 2 .
3.2.2. Active microwave data melt onset and total melt days.
The presence of liquid water greatly reduces snowpack volume scattering and increases microwave absorption, leading to a sharp reduction in active microwave backscatter coefficient sigma naught (σ 0 ) (Stiles and Ulaby 1980). This opposite response makes active microwave data an independent measure from the passive microwave for the remote NovZ and SevZ ice caps, where long term in-situ meteorological and glaciological data are limited for recent decades. Snowpack melt is detected when σ vv 0 drops below a commonly used threshold M, which is defined as M = σ wn 0 − ω where σ wn 0 the winter average of σ vv 0 (usually from January to February) and ω is a user-defined and region-specific constant (Trusel et al 2012, Wismann 2000. Sharp and Wang (2009)  or ω 2 = 3.5 dB from σ wn 0 for consecutive one or three days, respectively, are classified as melt days. Although we used the lower resolution 'egg' product (4.45 × 4.45 km 2 ) to minimize noise and maintain a comparable grid size with other data sets, melt records retrieved using the same two ω thresholds are consistent with earlier published results (Sharp and Wang 2009 , see supplementary notes and table S2). Therefore, this algorithm is applied to the entire lower resolution ('egg') QSCAT dataset to compare to the passive microwave results in the present work.
We resampled the 8.9 × 8.9 km 2 T b to four 4.45 × 4.45 km 2 QSCAT sub-grids based on nearest neighbor method and calculated the frequency distribution of TMD difference between SSM/I and QSCAT during the period of sensor overlap (2000-2007) ( figure S4(b)). QSCAT detected significantly more melt days than SSM/I. This is probably due to the fact that QSCAT data are sensitive to subsurface liquid water even when the ice surface is frozen (Hall et al 2009, Steffen et al 2004. Incomplete refreezing of water in deep snow layers and the capacity of firn to store water impede the active backscatter to return immediately to its pre-melt value when the air temperature decreases below 0°C, which questions the reliability of TMD estimated from active sensors (Rotschky et al 2011, Wismann 2000. In view of this, TMD is conservatively estimated from MERSL-BYU passive microwave data (SSM/I and AMSR-E) for 1995-2011. In accordance with the earlier definition, the active sensor-derived ice cap-wide MOD is the first day when more than 90% of analyzed pixels' σ 0 values drop below the threshold M (that is, 90% of the pixels concurrently experience melt). Through careful examination of backscatter variability compared to SSM/I and AMSR-E T b , thresholds of ω = 1.7 dB and ω = 1.0 dB for MOD detection are proposed for ERS and ASCAT, respectively. The thresholds agree well with SSM/I and AMSR-E algorithms (figure S5) and allow MOD to be retrieved for 1992-1994 and 2012 on both islands, extending the MOD time series to 1992-2012.
3.2.3. Partial correlation analysis. The correlation coefficient between sea ice extent and total snowmelt days as interpretation of the relationship alone can be misleading because sea ice and glaciers can equally co-respond to regional temperature change in the Russian High Arctic. Therefore, partial correlation analysis is adopted here, which calculates the correlation between two variables with the effect of a set of controlling variables removed. In this present research, the archipelago averaged 850 hPa geopotential height temperatures from NCEP-NCAR Reanalysis-1 is set as the controlling variable.

Regional temperature variations
The summer (June-September) mean 850 hPa geopotential height temperature time series from NCEP-NCAR and ERA-Interim are highly correlated over both NovZ (r = 0.952, pvalue < 0.001) and SevZ (r = 0.875, p-value < 0.001). ERA-Interim generally estimates a higher 850 hPa level Environ. Res. Lett. 9 (2014) 045009 M Zhao et al temperature and a slightly more pronounced temperature increase than NCEP-NCAR (figures 4(a), (b)). The temperature change from NCEP-NCAR during 1995-2011 are 0.59°C/decade (p-value = 0.11) and 0.80°C/decade (pvalue < 0.05) for SevZ and NovZ, respectively, while the two rates based on ERA-Interim increase to 0.85°C/decade (pvalue = 0.056) and 0.88°C/decade (p-value = 0.065), respectively. Over NovZ, temperatures from NCEP-NCAR have slightly stronger correlation with Malye Karmakuly stationmeasured average summer temperature on southwestern NovZ (r = 0.847, p-value < 0.001) than ERA-Interim (r = 0.830, p-value < 0.001). Over SevZ, both of the two reanalysis datasets have weak correlations with Golomjannyj Island station-measured summer mean temperatures. Golomjannyj Island station is located far west of SevZ into the Kara Sea with only 7 m a.s.l, and is strongly influenced by the cold ocean environment due to melting sea ice in summer. Opel et al (2009) found no correspondence between the melt-layer amount and the Golomjannyj station summer surface air temperatures, which further questions the basis for referencing Golomjannyj station temperature in glacier melting studies on SevZ. Additionally, our TMD results are strongly correlated with NCEP-NCAR reanalysis summer temperatures on both NovZ and SevZ. Therefore, we have confidence that the NCEP-NCAR reanalysis temperatures used in the present research are reliable and tend to make a more conservative estimate of the large-scale temperature change at 850 hPa geopotential height in NovZ and SevZ relative to ERA-Interim. Our decadal results confirm the strong positive relationship between mean TMD and the average June-September NCEP-NCAR reanalysis air temperature at geopotential height of 850 hPa, suggested by Sharp and Wang (2009) using a five year satellite record. The slopes of the linear regressions are ∼8 days°C −1 (r = 0.843, p-value < 0.0001, figure S6(a)) on SevZ and ∼10 days°C −1 (r = 0.787, pvalue < 0.0002, figure S7(a)) on NovZ from 1995 to 2011.

TMD correlation with late summer sea ice extent
Simple regression suggests that SevZ's TMD is significantly anti-correlated to Laptev Sea (r = −0.735, p-value < 0.001, figure S6(b)) and Kara Sea (r = 0.678, p-value < 0.003, figure  S6(d)) September sea ice extent, and NovZ's TMD is also significantly negatively correlated with annual late summer sea ice extent in the Barents Sea (r = 0.592, p-value < 0.02, figure S7 (b)) as well as that in the Kara Sea (r = 0.656, pvalue < 0.005, figure S7(d)). By setting the local temperature as the controlling variable, partial correlation result reveals that SevZ's TMD is still statistically anti-correlated to both Laptev Sea and Kara Sea late summer sea ice extent (table 1, figure S6(c), (e)). NovZ's TMD is also negatively related to Kara Sea September sea ice extent; however, it is not statistically related to Barents Sea late summer sea ice extent (table 1, figure S7(c), (e)).

Uncertainties
There are several uncertainties to consider in this study. Since we narrowed down our area of interest to ice cap interiors in order to avoid a mixed pixel effect, we maintained pixels that typically have later and shorter melt compared to lower elevation marginal ice (that may have been partially excluded), and the overall MOD and TMD are likely to be earlier and longer than the results presented here. Additionally, the SIR technique might smooth out some transient melt events when calculating the spatially weighed average of microwave signals (Long and Daum 1998). Therefore, the trends could be even more significant than reported here for each glacierized area. The switch between synthetic aperture radar and scatterometer mode of Advanced Microwave Instrument (AMI) onboard the ERS-1 and ERS-2 satellites resulted in many temporal coverage gaps. For the highest possible spatial resolution, 'all-pass' images are used by MERSL-BYU with a decrease in the temporal resolution to once per 5-6 days. Because of the low temporal resolution, it captured MOD with an uncertainty of 5-6 days. For instance, melt onset detected from ERS at the AWS station in 1999 was in the window of day 181-186 ( figure 2(a)). We set the MOD to be  June-September mean reanalysis temperature as the controlling variable. e-g Trend is determined by the slope of simple regression. a-g '**' and '*' indicate that the coefficient or trend is 95% and 90% statistically significant based on two-tailed Student's t-test, respectively.
the first day of the window, which might bias the MOD to a slightly earlier date. Because ERS is in the early part of our study period (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000), this approach makes the resulting trends more conservative (i.e. this does not artificially increase the trend toward earlier melt). In addition, MERSL-BYU ERS-1/2 data have a data gap from April 27th to June 3rd, 1996, and missed a significant melt event over NovZ on May 28th, 1996, which was detected by SSM/I data ( figure 3(b)). Therefore, we excluded ERS-derived melt onset for NovZ in 1996 from our analysis. Furthermore, MERSL-BYU AMSR-E data had 8 missing observations during the intense melt season in the year 2010. Those missing days were in the middle of melt and refreeze cycles and the contemporaneous ASCAT signals did not show an increase in σ 0 values on either NovZ or SevZ; therefore the present work assumes that melting occurred on those missing days and added 8 more days to each pixel's AMSR-E-derived 2010 TMD for both islands. In summary, we do not think these uncertainties jeopardize our conclusions.

Temperature control on snowmelt
The larger slope of the regression between TMD and temperature on NovZ indicates that NovZ is more sensitive to temperature variations at the 850 hPa level compared to SevZ. Sharp and Wang (2009) also discovered a decreasing sensitivity from NovZ to SevZ. Based on NCEP-NCAR reanalysis data, NovZ summer mean temperature at 850 hPa geopotential height increased at a significant rate of 0.80°C decade −1 (p-value < 0.05) during 1995 to 2011, while SevZ increased at a lower and less significant rate of 0.59°C y −1 (p-value = 0.11). The significant temperature increase and higher sensitivity to temperature change probably contributes to the larger ice mass loss rate (∼0.32 Mt a −1 km −2 ) of NovZ compared to SevZ (∼0.08 Mt a −1 km −2 ) between October 2003 and October 2009 (Moholdt et al 2012). The anomalously short melt season in 2010 decreases the overall statistical significance of TMD increase on NovZ; however, it is important to note that NovZ has experienced quite a few abnormally long melt years (2001, 2005, 2007, and 2011) in the past decade (figure 3 (c)). The interquartile range of TMD on NovZ increases at a rate of 0.43 d yr −1 (p-value < 0.01) during 1995-2011 (figure 3 (c)), indicating that variations in NovZ TMD are increasing. During 1995-2011, NovZ melted ∼23 more days than SevZ every year, which is consistent with the differing climate: SevZ is colder and drier than NovZ. The inherently longer melt season and the increasing snowmelt variations might also contribute to the recent higher ice loss on NovZ. Ice core data from Akademii Nauk ice cap on SevZ covering the period 1883-1998 indicate pronounced 20thcentury temperature changes in the Russian High Arctic, a strong rise in the early 1900s with the absolute temperature maximum in the 1930s (Opel et al 2009(Opel et al , 2013. The early 1900s warming is also documented in coastal marine sedimentary record on NovZ along eastern Barents Sea (Polyak et al 2004). Since temperature is a dominant factor that governs glacier summer melting, these pieces of evidence suggest a strong melting period during the early 1900s in Russian High Arctic. Further evidence for the existence of an intense melt period is the increase of melt-layer content in the Akademii Nauk ice core during the early 1900s (Opel et al 2009). The amount of melt layers on Akademii Nauk remains high until about 1970 and declines significantly since then, indicating a moderate melt period during 1970-1998. This pattern is consistent with a decrease of NCEP-NCAR reanalysis summer mean temperatures from 1960s to 1990s ( figure 4(c)). Meanwhile, our results show a significant increase in total snowmelt days on SevZ since 1996, coinciding with a strong recent decade temperature increase (figure 4(c)). Therefore, year 1996 might be a turning point from a cold period to a warm period for SevZ. However, reanalysis temperature data indicate large variations in glacier melting on NovZ and a pronounced recent decade of warming (figures 4(b), (c)).

Sea ice influence
Co-response to the regional temperature pattern is a straightforward interpretation of the negative relationship between late summer sea ice extent and glacier total snowmelt days in Russian High Arctic (e.g. Serreze et al 1995). After removing the large-scale temperature effect, partial correlation analysis suggests that glacier melt on SevZ is still statistically anti-related to Laptev Sea as well as Kara Sea sea ice extent, and NovZ glacier melt nevertheless significantly negatively correlate with Kara Sea sea ice extent (table 1). Sea ice decline in the Barents Sea seems to have limited influence on NovZ glacier melt. Reduced offshore sea ice extent, i.e. greater open-water fraction, can enhance onshore advection of sensible and latent heat fluxes (Rennermalm et al 2009). The insignificant late summer sea ice decline in Barents Sea during 1995 to 2011 might not affect the advection pattern and amount of heat flux from Barents Sea to NovZ (table 1), resulting in the less significant anti-relationship between NovZ TMD and Barents Sea sea ice extent.
Sea ice reduction can expose increasing areas of open water in summer to evaporation and change the large-scale atmospheric circulation, and has been shown to intensify summer precipitation over the Arctic and Antarctic as well as Northern Europe (Bamber et al 2004, Francis 2013, Screen 2013, Serreze et al 2012, Tietäväinen and Vihma 2008, Weatherly 2004. Sea ice influence on precipitation is also observed in Russian High Arctic. Opel et al (2009) analyzed stable water-isotope ratios from the uppermost 57 m of the ice core on the Akademii Nauk ice cap and found that sea ice cover is the main controlling factor for summer and autumn evaporation in the Kara Sea, and a lower sea ice extent allows higher evaporation rates and leads to larger contribution of regional moisture to the Akademii Nauk ice cap precipitation. Opel et al (2013) also discovered from Akademii Nauk sodium record that less sea ice cover and more open water in warmer times could increase the proportion of highly salted regional moisture, resulting in a more efficient transport and deposition of sea salt onto the Akademii Nauk ice cap.  Moholdt et al (2012) determined recent precipitation anomalies over NovZ and SevZ from the average and standard deviation of three global precipitation products: the Global Precipitation Climatology Project Version 2.2, the NCEP Climate Forecast System Reanalysis, and the ERA-Interim reanalysis. They found a slightly higher precipitation rate in 2004-2009 relative to the 1980-2009 mean, especially for Novaya Zemlya (89 ± 71 kg m −2 a −1 ), probably due to recent rapid sea ice reduction and increasing temperatures (Bamber et al 2004, Opel et al 2009, Serreze et al 2012. Although reanalysis precipitation products are likely to contain errors, the suggested positive precipitation anomaly is consistent with the spatial pattern of ice cap elevation change over NovZ in recent decades, whose ice cap interior has grown significantly while the periphery has thinned (Moholdt et al 2012). This pattern resembles the anomalous growth of central Austfonna ice cap in Svalbard during 1996-2002, where an increase in solid precipitation over ice cap center due to declining perennial sea ice extent in Barents Sea is the most likely explanation (Bamber et al 2004). Recent studies showed that NovZ and especially SevZ have experienced only slightly negative mass balance (Bassford et al 2006, Dowdeswell et al 1997, Gardner et al 2013, Matsuo and Heki 2013, Moholdt et al 2012, Zeeberg and Forman 2001; however, our study suggests that snowmelt becomes longer on both islands in the past decades, particularly significant on SevZ. The recent positive precipitation anomaly is very likely to be real and offset the increasing melt trend over both islands. In summary, the above discussion suggests that sea ice loss might enhance precipitation on NovZ and SevZ by exposing larger open ocean water to evaporation or amplifying regional temperature through ice-albedo feedback and water vapor greenhouse effect (

Conclusions
Our analysis determines glacier melt onset dates and total melt days on Novaya Zemlya and Severnaya Zemlya using multiple active and passive microwave sensors for the past two decades with a nominal grid size of 8.9 × 8.9 km 2 . Melt onset date has a strong earlier trend in Severnaya Zemlya and a highly variable recent record in Novaya Zemlya. Total melt days in both regions show an increase over time and are highly consistent with local temperature fluctuations. We discover a statistically significant negative relationship between regional late summer sea ice extent and glacier snowmelt, and discuss its possible influence on ice mass balance. This provides a potential explanation for the slightly negative mass balance rates on Novaya Zemlya and Severnaya Zemlya in the past several decades (Bassford et al 2006, Dowdeswell et al 1997, Gardner et al 2013, IASC 2010, Jacob et al 2012, Matsuo and Heki 2013, Moholdt et al 2012, Zeeberg and Forman 2001 despite recent warming in the Russian High Arctic (figure 4, IASC 2010). Our analysis suggests that the covariability pattern between sea ice retreat and a longer glacier snowmelt season is not simply a coresponse to regional temperature change. In combination with research on Greenland (Rennermalm et al 2009) and Svalbard (Bamber et al 2004, Rotschky et al 2011, shrinking sea ice is likely enhancing terrestrial ice snowmelt in a pan-Arctic perspective.