Introduction
The Antarctic Peninsula (AP) is a region that is warming far faster than the global average. During the second half of the 20th century the annual mean air temperature increased by 2.5–3.0 °C (Vaughan et al. Reference Vaughan, Marshall, Connolley, Parkinson, Mulvaney, Hodgson, King, Pudsey and Turner2003, Turner et al. Reference Turner, Colwell, Marshall, Lachlan-Cope, Carleton, Jones, Lagun, Reid and Iagovkina2005, Trenberth et al. Reference Trenberth, Jones, Ambenje, Bojariu, Easterling, Klein Tank and Parker2007, Steig et al. Reference Steig, Schneider, Rutherford, Mann, Comiso and Shindell2009). The number of days with positive air temperature and consequently the duration of melt events every year have increased at a rate of 0.5±0.3 days yr-1 over the 1980–2002 period (Torinesi et al. Reference Torinesi, Fily and Genthon2003). This climatic scenario increases surface melt, which accelerated over the 20th century (Abram et al. Reference Abram, Mulvaney, Wolff, Triest, Kipfstuhl, Trusel, Vimeux, Fleet and Arrowsmith2013). Since the millennium, a significant cooling trend was observed (Turner et al. Reference Turner, Lu, White, King, Phillips, Hosking, Bracegirdle, Marshall, Mulvaney and Deb2016) that also has glaciological implications (Oliva et al. Reference Oliva, Navarro, Hrbáček, Hernández, Nývlt, Pereira, Ruiz-Fernández and Trigo2017).
Surface melt has been regarded as one of the major factors associated with the breakup and disintegration of the Larsen A and B ice shelves (Scambos et al. Reference Scambos, Hulbe and Fahnestock2003). Surface runoff (SR) from land-based ice masses directly and indirectly influences short-term sea level changes. It adds freshwater to the oceans and lubricates the glacier bed, causing potential changes in, and acceleration of, the glacier mass flux (Pfeffer et al. Reference Pfeffer, Meier and Illangasekare1991, Zwally et al. Reference Zwally, Abdalati, Herring, Larson, Saba and Steffen2002, Vaughan Reference Vaughan2006, Osmanoglu et al. Reference Osmanoglu, Navarro, Hock, Braun and Corcuera2014). Recent studies have shown that ocean temperature is the leading factor that determines large ice shelves' thinning (Pritchard et al. Reference Pritchard, Ligtenberg, Fricker, Vaughan, van den Broeke and Padman2012). The resulting decreased buttress (Fürst et al. Reference Fürst, Durand, Gillet-Chaulet, Tavard, Rankl, Braun and Gagliardini2016), in turn, accelerates the ice flux of marine-terminating glaciers, resulting in a dynamic thinning much larger than the losses caused by surface melt (Wouters et al. Reference Wouters, Helm, Flament, Wessem, Van, Ligtenberg and Bamber2015). Though not the main contributors to mass loss and ice shelf breakup, surface melt and runoff are strongly associated with the observed changes on the AP.
Surface melt, surface runoff and ice discharge from glaciers by calving are the main drivers for freshwater input in the surrounding oceans of the AP. Freshwater supplies are important components in coastal ecosystems because they influence the physical and chemical settings of the water column, affecting both the structure and the function of coastal food webs (Moline et al. Reference Moline, Claustre, Frazer, Schofield and Vernet2004). Studies in shallow, circulation-restricted bays of King George Island showed that suspended particulate matter transported by glacier runoff can reduce light availability in shallow waters (Schloss et al. Reference Schloss, Abele, Moreau, Demers, Bers, González and Ferreyra2012) and governs the availability of macronutrients for phytoplankton production (Nędzarek Reference Nędzarek2008, Kienteca Lange et al. Reference Kienteca Lange, Tenenbaum, Tavano, Paranhos and Campos2014). Even far from the coast, Dierssen et al. (Reference Dierssen, Smith and Vernet2002) found phytoplanktonic blooms to be associated with glacier melt. Hodson et al. (Reference Hodson, Nowak, Sabacka, Jungblut, Navarro, Pearce, Ávila-Jimenez, Convey and Vieira2017) reported the occurrence of biologic activity hotspots associated with iron enrichment supplied by glacial weathering.
Previous studies investigated surface melt and runoff on AP primarily through remote sensing and regional climate modelling. Using passive microwave data (1978–present), studies addressing changes in melt extent, duration and water quantities could not find a significant temporal trend. However, these studies note decreasing, though statistically insignificant, trends for melt extent, melt index and total annual melt (Torinesi et al. Reference Torinesi, Fily and Genthon2003, Liu et al. Reference Liu, Wang and Jezek2006, Kuipers Munneke et al. Reference Kuipers Munneke, Picard, Van Den Broeke, Lenaerts and Van Meijgaard2012, Trusel et al. Reference Trusel, Frey, Das, Munneke and Van Den Broeke2013). The only positive trend in melt duration was found by Torinesi et al. (Reference Torinesi, Fily and Genthon2003) from 1980–2000. They found that melting lasts for an average of 50 days on the AP, reaching a maximum of 100, with an increasing trend of 0.5±0.3 day yr-1 (1980–1999). Liu et al. (Reference Liu, Wang and Jezek2006) determined a median melt duration of 59 days with an absolute variation of 5.39 days and a melt extent ranging from 2.5x105–3.5x105km2. Barrand et al. (Reference Barrand, Vaughan, Steiner, Tedesco, Kuipers Munneke, Van Den Broeke and Hosking2013) combined remote sensing observations and simulations from the regional model RACMO2 to investigate melt conditions, finding melt extents ranging from 2.8x105–3.3x105 km2 (QuickSCAT) and 1.8x105–3 x105 km2 (RACMO2 regional atmospheric model) from 2000–2009.
Quantitative estimates of melt rates or freshwater input to the ocean are even scarcer. Van de Berg et al. (Reference Van de Berg, van den Broeke, Reijmer and van Meijgaard2005) showed a maximum annual melt of 0.5 m water equivalent (w.e.) yr-1 occurring on the Larsen Ice Shelf's northern edge from 1958–2002. Vaughan (Reference Vaughan2006) estimated values for surface melt in the year 2000 of 54±26 Gt (450±216 mm w.e.). He considered only the grounded ice portion of the AP and predicted an increase in surface melt to 100±46 Gt by 2050, assuming the current increasing trend of annual mean air temperature. Hock et al. (Reference Hock, De Woul, Radic and Dyurgerov2009) estimated that from 1961–2004, an increase in the sea level of 0.22±0.16 mm yr-1 resulted from the contribution of mass loss from glaciers and ice caps surrounding Antarctica, mostly AP. More recently, Kuipers Munneke et al. (Reference Kuipers Munneke, Picard, Van Den Broeke, Lenaerts and Van Meijgaard2012) found surface melt values ranging from 20–130 Gt yr-1 for the AP. The analysis considered the entire Antarctic ice sheet and surrounding ice shelves, from 1979–2009, and showed no statistically significant decreasing trend.
This study presents comprehensive estimates for surface melt and runoff on the Antarctic Peninsula and surrounding islands. It uses a positive degree day-based glacier melt model adapted from Vaughan (Reference Vaughan2006) and Pfeffer et al. (Reference Pfeffer, Meier and Illangasekare1991), driven by data from the global reanalysis project ERA-Interim (ERAI) from the European Centre for Medium-range Weather Forecast (ECMWF). The model is calibrated and validated via a multi-criteria scheme using in situ observations and long-term surface mass balance (SMB) records. The model covers the period 1981–2014.
Study site
The AP is a region with strong latitudinal and longitudinal gradients of climatic parameters (Morris & Vaughan Reference Morris and Vaughan2003). The central mountain range, reaching more than 2000 m a.s.l., forms a major obstacle in the southern hemisphere polar vortex (Fig. 1). The climate of the AP west coast and adjacent islands is cold maritime, while the conditions are much more continental on the east coast and further south. Consequently, summer melt regularly occurs along the west coast at lower elevations. Frequent warm foehn-type winds on the east coast are known to cause considerable surface melt, and cold barrier winds lead to cold air mass outbreaks from the Filchner-Ronne Ice Shelf along the AP's mountains to the north. Such foehn events impact the surface temperature regionally in timescales varying from hours to seasons (Cape et al. Reference Cape, Vernet, Skvarca, Marinsek, Scambos and Domack2015). Larsen C&D, Wilkins and George VI ice shelves are the largest low-elevation areas. The plateau of the Antarctic Peninsula is dominated by a dry snow zone, with mean annual surface temperatures below -11°C (Rau & Braun Reference Rau and Braun2002). Large mountain glaciers drain from the plateau to both sides of the AP. The small islands surrounding the peninsula are covered by ice caps and ice fields.
Fig. 1 Map showing the study area and data-source distribution. Elevation contours are extracted from RAMP-DEM (Liu et al. Reference Liu, Jezek, Li and Zhao2015). Weather stations are represented by yellow diamonds and labelled with numbers, as shown in Table II. The top right detailed map of King George Island shows the glaciological stakes and ultrasonic ranger used for calibration and validation in the Bellingshausen Dome and the Fildes Peninsula. Magenta squares show locations in Livingston and Vega islands, where long-term surface mass balance records are available.
Data
The AP does not have a well-distributed weather station network with long time series (see Fig. 1). Under these circumstances, a feasible alternative to the spatial interpolation of uneven observations is to use data provided by reanalysis projects. Such datasets incorporate available observations to minimize the errors of the prediction model.
ERAI is a reanalysis product from ECMWF. It provides data from 1979–present with a finest spatial resolution of approximately 0.25ox0.25o. Compared with previous versions, such as ERA-40 and ERA-15, the temperature bias over Antarctica has been reduced. The orography in the model is an average of finer resolution digital elevation models, such as GTOPO30 (Gesch et al. Reference Gesch, Verdin and Greenlee1999). This averaging smooths the highly complex topography of the AP.
To better represent the topography, the RAMP DEM (Liu et al. Reference Liu, Jezek, Li and Zhao2015) was used as the input for an altitudinal lapse-rate downscaling method, which was evaluated using near-surface air temperature data from 28 weather stations. The resulting grid cell size after downscaling was 200 m x 200 m.
Direct ablation measurements were available from three sources: 1) Measurements in 10–14 day intervals at 29 mass balance stakes during the summer from 2007–2012 at the Bellingshausen Dome in King George Island, South Shetland Islands (stake line 1) (Table I; Mavlyulov Reference Mavlyulov2014), 2) continuous surface melt measurements available from a sonic ranging sensor (SR50) operated over 42 days of the 1997/98 summer at daily intervals (Braun et al. Reference Braun, Saurer, Vogt, Simoes and Gossmann2001, Braun & Hock Reference Braun and Hock2004), 3) mass balance stakes readings during summer field campaigns (stake line 2). These measurements were conducted during two summer seasons (1997/1998 and 1999/2000) on the Bellingshausen Dome (Braun et al. Reference Braun, Saurer, Vogt, Simoes and Gossmann2001, Braun & Hock Reference Braun and Hock2004).
Table I Data used for the calibration and validation of the surface melt and runoff model.
KGI denotes King George Island, LIV: Livingston Island, VI: Vega Island, Cal: calibration data, Val: validation data.
Spatially-integrated summer and winter mass balances of the Hurd and Johnsons glaciers in Hurd Peninsula, Livingston Island were available for ten years (Navarro et al. Reference Navarro, Jonsell, Corcuera and Martín-Español2013, data provided by the World Glacier Monitoring System). A 14-year record of the annual surface mass balance of Glaciar Bahía del Diablo, Vega Island was used as a quantitative reference for the north-eastern AP region (Skvarca et al. Reference Skvarca, De Angelis and Ermolin2004, Marinsek & Ermolin Reference Marinsek and Ermolin2015).
Methods
ERA-Interim temperature downscaling
To account for the spatiotemporal variation in the air temperature altitudinal lapse rates along the study area, the temperature and geopotential height differences between the 1000 hPa and 750 hPa pressure levels were computed. Thereafter, the ratio between the two differences for each day and grid element was calculated, and the given lapse rate (Eq. (1)) was applied to the elevation difference between the RAMP-DEM and the ERAI geopotential height at the surface (Eq. (2)). The result was then reduced from the temperature at 2 m (Eq. (3)).
$$lr\left( {\varphi ,\lambda ,d} \right){\equals}{\rm }{{T_{{1000}} \left( {\varphi ,\lambda ,d} \right){\minus}T_{{750}} \left( {\varphi ,\lambda ,d} \right){\rm }} \over {h_{{1000}} \left( {\varphi ,\lambda ,d} \right){\minus}h_{{750}} \left( {\varphi ,\lambda ,d} \right){\rm }}}$$
$$dh{\equals}h_{{ERA{\minus}sfc}} {\minus}h_{{RAMP{\minus}DEM}} $$
and
$$\equaligno\hskip-1pt { T_{{DS}} \left( {\varphi ,\lambda ,d} \right){\equals}{\rm }T_{{2m}} \left( {\varphi ,\lambda ,d} \right){\plus}(dh\left( {\varphi ,\lambda ,d} \right)\,{\rm {\asterisk}}\,lr\left( {\varphi ,\lambda ,d} \right)),\vskip5pt $$
where lr is the lapse rate, T 1000, T 750, h 1000 and h 750 are the temperatures and geopotential heights at 1000 hPa and 750 hPa, respectively, h ERA-sfc is the geopotential height of the surface, h RAMP-DEM is the elevation given by RAMP DEM (Liu et al. Reference Liu, Jezek, Li and Zhao2015), T DS is the temperature downscaled, T 2m is the temperature 2 m above the surface, φ and λ are the latitude and longitude, and d is the day. For each model grid cell, the ERAI cell with the cell centre closest to the model grid centre was used.
It is worth noting that the lapse-rate method for air temperature downscaling only accounts for the decrease of the air temperature due to the vertical variation inside the ERA-Interim grid elements. It does not include any further dynamic processes that may affect the air temperature which are not already represented in the ERA-interim global climatic model. It also represents the free atmosphere and does not account for any local effects or boundary-layer processes.
Surface melt and runoff computations
The positive degree day (PDD) is the sum of daily positive near-surface air temperatures during a given period (in this case, one year) (Vaughan Reference Vaughan2006). When multiplied by a melt factor (m), it provides the total surface melt (M):
$$M\left( {\phi ,\lambda ,\tau } \right){\equals}m \mathop \sum\nolimits_{d{\equals}1}^{d{\equals}365} T_{{DS}} \left( {\phi ,\lambda ,d} \right)\alpha \left( {\phi ,\lambda ,d} \right)$$
$$\left\{ {\matrix{ {\alpha \left( {\phi ,\lambda ,d} \right){\equals}1,if{\rm }T_{{DS}} \left( {\phi ,\lambda ,d} \right)\,\gt\,0} \cr {\alpha \left( {\phi ,\lambda ,d} \right){\equals}0,if{\rm }T_{{DS}} \left( {\phi ,\lambda ,d} \right)\,≤\,0} \cr } } \right.$$
where d=1 is the first day of the annual melting periods (assumed to 1 October in this case) and τ is the year.
The energy balance at the surface implies that melt can occur when the air temperature is equal to or below 0 °C, and even not occur when the temperature is positive. In this sense, the 0 °C threshold for melt occurrence may be considered an artificial threshold. In order to determine the more suitable threshold, the model was used to consider five different thresholds for the same period covered by the work of Barrand et al. (Reference Barrand, Vaughan, Steiner, Tedesco, Kuipers Munneke, Van Den Broeke and Hosking2013), i.e. from 1999−2002. The mean melt duration of the period as estimated by our model was then compared against the QuickSCAT-derived estimates obtained by Barrand et al. (Reference Barrand, Vaughan, Steiner, Tedesco, Kuipers Munneke, Van Den Broeke and Hosking2013). The classical PDD approach considers that melt occurs when the air temperature is positive. The work presented here also explores melt when the air temperature is equal to and above 0 °C, above -1°C, above -0.5 °C, above 0.5 °C and above 1 °C. The latter two thresholds were tested in order to account for the known overestimation of the ERAI-derived air temperature on the AP.
Surface runoff (Q) is calculated as the difference between surface melt and the amount of meltwater retained in the snowpack through refreezing, pore filling and capillarity (M 0):
$$Q{\equals}M{\minus}M_{0} $$
with M 0 given by (Pfeffer et al. Reference Pfeffer, Meier and Illangasekare1991) as
$$M_0 {\equals}{c \over L}CT_{f} {\plus}\left( {C{\minus}M} \right)\left( {{{\rho _{{pc}} {\minus}\rho _{c} } \over {\rho _{c} }}} \right),$$
where c is the heat capacity of ice, L is the latent heat of fusion of ice, C is the snow accumulation, T f is the temperature of the firn at the beginning of the melt season (in positive Celsius degrees below the freezing point), ρ c is the initial firn density (taken as 400 kg m-3) and ρ pc is the pore close-off density (taken as 830 kg m-3). After setting ρ c =400 kg m-3 and ρ pc =830 kg m-3 (Vaughan Reference Vaughan2006), one gets:
$$M_0 {\equals}\left( {0.003T_{f} {\plus}0.52} \right)C.$$
After assuming that T f is the mean air temperature of the previous year in each grid cell, the only unknown variable for calculating M 0 is C. C is approximated as the annual accumulation for each grid element given by ERAI (resampled by nearest neighbour to the model grid). Hence, the accumulation and temperature-driven changes in retention for each year is estimated, and the surface runoff is estimated using Eq. (5).
Calibration and validation of the model
A multi-criteria calibration and validation scheme was applied, with rigorous data splitting for model calibration, data downscaling and model validation.
To validate the downscaling approach, monthly PDDs that were computed from the downscaled ERAI data and from weather station data (records longer than ten years) were compared. The analysis was restricted to the months from October to March to focus on the melt periods. The PDDs' correlation strength and the sum of errors were utilized as quality indicators.
The melt model was calibrated by tuning the melt factor using ERAI downscaled data in Eq. (4) to achieve an optimal fit to the ablation records (N=300) of stake line 1. A melt factor of 5.4 mm w.e. K d-1 was determined with r 2=0.65 using all available records. When stakes are considered individually, the melt factor takes values between 2.2 and 12.6 mm w.e. K d-1 (although both of these are extreme values). The melt factor is well within the range reported for comparable Arctic regions and at the upper bound for Antarctica (Huybrechts & Oerlemans Reference Huybrechts and Oerlemans1990, Braithwaite & Zhang Reference Braithwaite and Zhang2000, Hock Reference Hock2003).
Independent validation of the model performance was carried out using a variety of datasets at different locations and different time intervals. The melt rates at stake line 2 and the SR50 continuous record were compared over a six-week melt period, both on King George Island. The surface mass balance records provided an integrated error estimate for summer and winter (Livingston Island) (Navarro et al. Reference Navarro, Jonsell, Corcuera and Martín-Español2013) or an entire glacier mass balance year (Vega Island) (Skvarca et al. Reference Skvarca, De Angelis and Ermolin2004, Marinsek & Ermolin Reference Marinsek and Ermolin2015).
Results
PDD threshold and spatial representation of the mean melt duration (MMD)
By spatially comparing the annual mean melt duration from 1999–2002 (Fig. 2), it can be observed that the best performance is achieved when using the 0 °C threshold. The energy balance at the surface determines that melting can occur at temperatures below 0 °C. However, the comparisons made between the ERAI-derived and the weather stations-derived PDDs (see below and Table II) show that both are very well correlated, although the absolute values are overestimated. This may explain the best performance of the 0 °C threshold for determining the mean melt duration, although it could be expected that a lower temperature would result in a better representation of the mean melt duration.
Fig. 2 Comparison of the mean melt duration (MMD) derived from ERAI-based air temperature data against QuickSCAT-derived maps. Six different thresholds are applied to ERAI-derived data to determine the occurrence of melt: a. T>1°C, b. T>0.5°C, c. T>0°C, d. T≥0°C, e. T>-0.5°C, and f. T>-1°C.
Table II Weather stations used in this study for the validation of PDDs estimated from ERAI. The correlation coefficient and RMSE were calculated between monthly PDDs derived from weather station (WS) measurements and ERAI estimation, excluding the months from April to September.
Downscaling and melt model performance
The ERAID-derived PDDs showed very good agreement with weather station data. Correlation coefficients ranged from 0.58–1.0, with the majority above 0.9 (Table II). Nevertheless, absolute values were overestimated by ERAID-derived data.
The surface melt was validated using two other independent point ablation measurements carried out in the Bellingshausen Dome (stake line 2 and SR50). Compared to these datasets, the modelled surface melt produced a slight underestimation of the surface ablation. After 41 days of ablation monitoring with the SR50, a final difference of 80 mm w.e. was found between the modelled and the measured ablation. The root mean square error (RMSE) was 7.19 mm w.e., and the bias was 4.9 mm w.e. The linear correlation of the ablation estimates with stake line 2 was 0.86 (Fig. 3), with RMSE of 19.08 mm w.e. and bias of 7.14 mm w.e.
Fig. 3 a. MF estimation, b. validation against ablation measured at stake line 2 (1997/1998, 1999/2000) and c. SR50 (40 consecutive days, from December 1997 to January 1998).
The surface melt estimations in the Hurd Peninsula (Livingston Island) were compared to measurements carried out by Navarro et al. (Reference Navarro, Jonsell, Corcuera and Martín-Español2013) on the Hurd Glacier (HG) and Johnsons Glacier (JG) from 2002–2011 (Fig. 4). The same study also measured the winter accumulation, which was used to evaluate the accumulation estimate over the same periods. The correlation of the modelled ablation was 0.61 with HG and 0.68 with JG, with RMSE of 158 mm w.e. The accumulation correlations were 0.8 (HG) and 0.49 (JG), with RMSE of 247 mm w.e. Values always ranged on the same scale and did not show a fixed over- or underestimation. Additionally, the annual surface mass balance of Glaciar Bahía del Diablo (Vega Island) was estimated and compared to field measured data from this site (see Fig. 4). The error was lower than 200 mm w.e. for nine of the 13 analysed years, the correlation between the measured and the modelled SMB was 0.67, and the RMSE was 214 mm w.e.
Fig. 4 Comparison of modelled area-averaged SMB in: a. the Hurd Peninsula (Livingston Island) with summer and winter SMB from the Hurd (HG) and Johnsons (JG) Glaciers (Navarro et al., Reference Navarro, Jonsell, Corcuera and Martín-Español2013) and b. Bahía del Diablo Glacier's (Vega Island) annual SMB (Skvarca et al., Reference Skvarca, De Angelis and Ermolin2004; Marinsek and Ermolin, Reference Marinsek and Ermolin2015).
The average ablation was underestimated in the SR50 and stake line 2 (SL2) validation site, by respectively 0.14 mm w.e. (2% below the observation) and 11.4 mm w.e. (20% above the observation). The average summer ablation of Livingston Island was overestimated by 120.45 mm w.e. (17% above the observation), whereas the winter accumulation was overestimated by 107.69 mm w.e. (16% above the observation). The annual mass balance of the Bahía del Diablo glacier was overestimated by 30.52 mm w.e. (20% above the observation). The best correlation between the modelled and the observed time series was found with the SL2 observations, which is expected because the SL2 is located at the melt factor calibration site. Although the SR50 measurements were also performed at the calibration site, the PDD method is known for its poorer performance on the daily timescale when compared to longer timescales.
The observational data available for validating this study consists of measurements performed in time intervals ranging from daily to annual, implying different magnitudes of ablation, accumulation and mass balance. To account for the diverse measurement magnitudes, the normalized RMSE was computed while considering the range of the observed measurements. The normalized RMSE ranged from 0.16–0.25, with the lowest value associated with the SL2 validation site and the highest value associated with the winter accumulation on Livingston Island. It was expected that the model would have a worse performance away from the melt factor calibration site because melt factors typically show a large variation in regions as large as the Antarctic Peninsula. However, a better comparison of errors would require longer measurements carried out at similar timescales.
Surface melt and runoff time series
The annual time series for the mean near-surface air temperature and PDDs and the total surface melt (M), M 0 and surface runoff (Q) are presented in Fig. 5 for the total, grounded and floating areas of the AP. The years correspond to the melt season that started in October of the previous year.
Fig. 5 a. Mean temperature, b. mean PDDs, c. total surface melt (M), d. total potential retention (M 0) and the total surface runoff (Q) annual time series for the total (continuous line, diamonds), grounded ice (dotted line, upward-pointing triangles) and floating ice (dotted line, downward-pointing triangles) areas. In d, the right axis shows scale for the floating area.
The mean near-surface air temperature does not represent the PDDs' interannual variability. Both series show r 2=0.11, which is not surprising considering that there are no PDDs during most months that compose the annual mean temperature.
Total surface melt had maximum values in 1985 (129 Gt) and 1993 (127 Gt) and a minimum value in 2014 (26 Gt). The mean and standard deviation values for the entire period and area were 75±54 Gt, or 46±15 Gt and 25±8 Gt when separately considering floating and grounded areas, respectively. Surface runoff had a maximum value in 1993 (40 Gt) and a minimum value in 2014 (0.37 Gt), with a mean of 9±8 Gt and, for floating and grounded areas, 5±6 Gt and 3±2 Gt, respectively.
Since 2008, both surface melt and runoff have persistently shown negative anomalies when considering the entire period. The same happened from 1996–2002 but only for surface runoff. The period 1996–2002 is also characterized by a lower interannual variability in the mean air temperature, PDDs and surface melt, whereas M 0 has not shown the same behaviour. Hence, the decoupling of surface runoff and surface melt for this period might be associated with high accumulation rates.
The spatial distribution of the mean and standard deviation of surface melt and runoff is presented in Fig. 6. The effects of foehn winds are clearly reflected in the mean surface melt (Fig. 6a) for the western border of the Larsen C Ice Shelf, showing that these types of events are large and persistent enough to be represented in a global reanalysis.
Fig. 6 Maps showing: a. mean surface melt (
$\bar{M}$
) and its standard deviation (σ
M
, a') over the 1980–2014 period and b. the mean surface runoff (
$\bar{Q}$
) and its standard deviation (σ
Q
, b') over the same period.
Surface melt and runoff regions
The study divided the AP in 11 sub-regions (Fig. 7) to analyse the spatial surface melt and runoff variations in the major drainage basins and ice shelves, bays and topographical features. Because the areas (A) of the regions considered are different, the specific surface melt and runoff (M/A and Q/A, in m w.e.) were computed. The largest specific surface melt and runoff correspond to the South Shetland Islands (SSI). Excluding the islands, specific surface melt and runoff are larger in the E-LB region.
Fig. 7 Specific surface melt (M/A) and surface runoff (Q/A) over the Antarctic Peninsula: a. to k. show the M/A (continuous line, circles) and Q/A (continuous line, squares) time series separated for each AP region; l. shows the mean and standard deviation of the surface runoff maximum elevation in each region; map m. delimits each region.
The regional time series have a higher correlation than the specific surface runoff time series. The specific surface melt and runoff of most areas are better correlated with their immediate neighbours to the east or the west sides. The northern areas (W-N1, W-N2, E-N and E-LB; Fig. 7) and northern islands (SSI and Joi) show very similar specific surface melt and runoff temporal behaviours; the differences strongly reflect the elevation profiles. E-N shows a higher correlation with SSI (0.91) and Joi (0.91) than W-N1 (0.86 and 0.74, respectively). E-LC correlates better with E-LB, although the absolute amounts are considerably lower for E-LC as compared to E-LB. E-S correlates better with W-S than any other region, showing that there is a more pronounced climatic difference, in terms of interannual variability, between the regions in the E-LC latitudinal range.
The distribution of the surface melt and surface runoff along the elevation profiles reveals that more than 90% of the specific surface melt (specific surface runoff) is produced at elevations between 0 m and 500 m (400 m). This proportion decreases with elevation in all regions, except for the E-LC area. In this area, SM/A is higher in the elevation range of 100–200 m than in 0–100 m. The occurrence of foehn-type winds in this area leads to higher temperatures near the border with the grounding line, increasing melt.
Specific surface runoff is better distributed along the profile, and 90% occurs below an altitude of 400 m. The SSI and the Joi also show similar distribution patterns. The W-Wi and the W-S have 45% and 50% of specific surface melt, respectively, restricted to 0–100 m. The W-MB shows 90% of specific surface runoff restricted to 0–200 m. The greatest difference is found between the northern regions (E-N, E-LB, W-N1, W-N2, SSI and Joi) and the southwestern regions (W-S, W-MB and W-Wi). E-S has a distribution with characteristics between these extremes. E-LC is a special case, where more melt is produced between 100–200 m a.s.l. than between 0–100 m a.s.l., as a result of foehn-type winds.
Discussion
The high correlation between ERAID-derived and weather stations-derived PDDs demonstrates the consistency of the temporal variability of the modelled PDDs. On the other hand, very high RMSEs show that the PDDs' absolute values are not well represented by ERAID. An overestimation of air temperature over Antarctica has been recognized in previous ERA datasets (van de Berg et al. Reference Van de Berg, van den Broeke, Reijmer and van Meijgaard2005), which is also observed in this study.
The comparative analysis between the ERAI-derived mean melt duration and the QuickSCAT-derived mean melt duration (Barrand et al., Reference Barrand, Vaughan, Steiner, Tedesco, Kuipers Munneke, Van Den Broeke and Hosking2013) showed that the first is underestimated on high-slope areas and overestimated on the northern tip of the Larsen C Ice Shelf. The terrain slope is a major difficulty for atmospheric models and even for remote sensing data processing. Therefore, a worse performance in high-slope areas is expected. Aside from the observation that the northern Larsen C mean melt duration overestimation occurs in flat areas, it may also be a topography-related problem. The north-western Larsen C is influenced by foehn-type winds, whose proper modelling is highly dependent on accurate representation of terrain elevation.
However, it is also worth noting that Luckman et al. (Reference Luckman, Elvidge, Jansen, Kulessa, Kuipers Munneke, King and Barrand2014) compared Envisat/ASAR melt duration maps with the same QuickSCAT-derived data (Barrand et al., Reference Barrand, Vaughan, Steiner, Tedesco, Kuipers Munneke, Van Den Broeke and Hosking2013), and found melt durations up to 25 days longer than the latter. The authors attribute this difference to the different spatial resolution of both datasets in the region of Larsen C influenced by foehn-type winds. This suggests that perhaps the QuickSCAT reference data is actually underestimating the melt duration in Larsen C, which would indicate that the mean melt duration found in the present study is less overestimated.
The approach used here for tuning the melt factor to ERAID PDDs using in situ ablation measurements allowed for a compensation of the PDDs' overestimation. Both the surface melt temporal variability and absolute surface melt values are consistent with the SR50 measurements. The 80 mm w.e. underestimation over 40 days is acceptable, particularly considering that a simple approach applied to global reanalysis data was used. A small underestimation and good correspondence with SL2 ablation measurements were also found. The comparisons with integrated SMB from Hurd Peninsula and Glaciar Bahía del Diablo show that the results presented here are reliable for a regional-scale analysis on the northern AP. However, it should also be noted that there was no quantitative melt or surface mass balance data to assess the quality of the results for more southern sites on the AP.
A remaining bias from the PDD overestimation not compensated for by the melt factor tuning would result in a larger melt area. However, the spatial distribution of PDDs is in agreement with recent work addressing the melt occurrence in AP through QuickSCAT and regional modelling (Kuipers Munneke et al. Reference Kuipers Munneke, Picard, Van Den Broeke, Lenaerts and Van Meijgaard2012, Barrand et al. Reference Barrand, Vaughan, Steiner, Tedesco, Kuipers Munneke, Van Den Broeke and Hosking2013, Trusel et al. Reference Trusel, Frey, Das, Munneke and Van Den Broeke2013, van Wessem et al. Reference van Wessem, Ligtenberg, Reijmer, van de Berg, van den Broeke, Barrand, Thomas, Turner, Wuite, Scambos and van Meijgaard2016). In some regions, the data in this study indicate a smaller melt area, which can be attributed to the enhanced resolution of the grid after the ERAI temperature downscaling. QuickSCAT data have spatial resolution of 5 km and RACMO2.3 of 5.5 km, whereas temperatures in this work were downscaled to a 200 m x 200 m grid.
Both surface melt and surface runoff spatial distribution and time series are generally in good agreement with the most recent studies addressing surface melt patterns in the AP (Tedesco & Monaghan Reference Tedesco and Monaghan2009, Kuipers Munneke et al. Reference Kuipers Munneke, Picard, Van Den Broeke, Lenaerts and Van Meijgaard2012, Barrand et al. Reference Barrand, Vaughan, Steiner, Tedesco, Kuipers Munneke, Van Den Broeke and Hosking2013, Trusel et al. Reference Trusel, Frey, Das, Munneke and Van Den Broeke2013, Välisuo et al. Reference Välisuo, Vihma and King2014, van Wessem et al. Reference van Wessem, Ligtenberg, Reijmer, van de Berg, van den Broeke, Barrand, Thomas, Turner, Wuite, Scambos and van Meijgaard2016), but the absolute values are higher than those presented by previous studies. Nevertheless, comparisons of surface melt derived from RACMO with 27 km resolution against QuickSCAT derived data showed an underestimation of SM, melt onset date and melt season duration by RACMO (Kuipers Munneke et al. Reference Kuipers Munneke, Picard, Van Den Broeke, Lenaerts and Van Meijgaard2012, Barrand et al. Reference Barrand, Vaughan, Steiner, Tedesco, Kuipers Munneke, Van Den Broeke and Hosking2013, Trusel et al. Reference Trusel, Frey, Das, Munneke and Van Den Broeke2013). Van Wessem et al. (Reference van Wessem, Ligtenberg, Reijmer, van de Berg, van den Broeke, Barrand, Thomas, Turner, Wuite, Scambos and van Meijgaard2016) present estimates of surface melt provided by RACMO with 5.5 km resolution that are considerably smaller than those obtained by the previous model version. Comparisons of QuickSCAT-derived melt area on Larsen C against higher resolution Envisat-ASAR derived data further suggest that even QuickSCAT-derived melt area is underrepresented in the region. Such discrepancies between the different studies depict the high uncertainties involved in modelling mass balance terms in the Antarctic Peninsula.
The effects of foehn-type winds are visible in the distributed mean surface melt, contrasting with their poor representation in previous modelling studies (Kuipers Munneke et al. Reference Kuipers Munneke, Picard, Van Den Broeke, Lenaerts and Van Meijgaard2012, van Wessem et al. Reference van Wessem, Ligtenberg, Reijmer, van de Berg, van den Broeke, Barrand, Thomas, Turner, Wuite, Scambos and van Meijgaard2016). In general, the surface melt time series presented here agrees well with the temporal behaviour of the surface melt time series estimated by Kuipers Munneke et al. (Reference Kuipers Munneke, Picard, Van Den Broeke, Lenaerts and Van Meijgaard2012) for the AP. In spite of the uncertainties associated with the approach used in the work discussed here, the results are consistent with direct field measurements (Figs 2 & 3).
Wouters et al. (Reference Wouters, Helm, Flament, Wessem, Van, Ligtenberg and Bamber2015) found a total mass loss of 300 Gt from 2000–2014. The study comprised the grounded portion of the AP southern sector. During the same period, and considering the grounded area of the entire AP, losses by surface runoff (SR) were 95 Gt. The contributions of the mass losses by SR to glacier thinning and possible sea level rise are thus small compared to ice flux acceleration and increased calving. Nevertheless, the role of the floating areas as a direct freshwater source to the adjacent ocean is remarkable. While they comprise only 24% of the total area, they produce 68% of the mean SR and 61% of the mean SM. Additionally, the fate of the liquid water is different from the fate of the ice when entering the ocean and it impacts the physico-chemical characteristics and the biota of the water differently.
The AP floating area has decreased by 18% since 1950 (Cook & Vaughan Reference Cook and Vaughan2010). This reduction is progressively advancing southwards. The sea level rise around the AP, predominantly driven by steric expansion (Rye et al. Reference Rye, Naveira Garabato, Holland, Meredith, George Nurser, Hughes, Coward and Webb2014), suggests that the oceanic forcing will continue to lead to further and more dramatic ice shelf disintegration. In the future, it is likely that the melting area of the AP will be further reduced. Both Larsen C and the Wilkins Ice Shelf are the main surface melt and surface runoff sources and are partially unstable. Jansen et al. (Reference Jansen, Luckman, Cook, Bevan, Kulessa, Hubbard and Holland2015) predicted that the Larsen C would undergo a large calving event caused by the development and propagation of a rift; this started on 12 January 2017 and reduced its area to approximately 90% of the pre-collapse size (Hogg & Gudmundsson Reference Hogg and Gudmundsson2017). After a large breakup event in February 2008, the 3100 km2 of the northern portion of the Wilkins Ice Shelf was at risk of collapse (Braun et al. Reference Braun, Humbert and Moll2009). A narrow ice bridge collapsed later, in April 2009 (Humbert et al. Reference Humbert, Gross, Müller, Braun, Van De Wal, Van Den Broeke, Vaughan and Van De Berg2010). Consequently, ice shelf breakups may reduce the absolute amounts of water input to the ocean in the form of surface runoff.
The accumulation rate on the AP has doubled since 1850 (Thomas et al. Reference Thomas, Marshall and McConnell2008), which, according to the Pfeffer et al. (Reference Pfeffer, Meier and Illangasekare1991) approach, increases meltwater retention. Nevertheless, Abram et al. (Reference Abram, Mulvaney, Wolff, Triest, Kipfstuhl, Trusel, Vimeux, Fleet and Arrowsmith2013) showed that surface melt intensified much faster. By their analysis, the melt intensity recorded in the ice core increased from 0.5% to 4.9% over the 20th century. The grounded area presently retains an average of 91% of the SM, which will reduce if the observed trends continue; this may lead to an increase in the maximum height where surface runoff occurs and an increase in SR intensity at lower elevations.
It is possible to differentiate the specific surface melt and runoff time series in the 11 regions, even though the model is fed by global reanalysis data. The temperature downscaling using altitudinal lapse-rate variations in space and time allowed for an East–West differentiation in the northern tip of the AP. The W-N2 region lies on the same latitudinal range of E-LB, and it shows specific surface melt and runoff approximately 50% lower than the latter region. The difference may be explained by topographical control (E-LB has a mean elevation 330 m lower) and smaller altitudinal temperature lapse rates on the east side of the AP (Morris and Vaughan, Reference Morris and Vaughan2003; this study).
The regional analysis of the altitudinal distribution of the specific surface melt and runoff revealed that they are determined, concomitantly, by the latitudinal range and East−West climatic differences. Though the western side is warmer, in the southwestern regions the specific surface melt and runoff are less uniformly distributed over the topographic profile due to the cold climate and higher temperature altitudinal lapse rates. Due to higher accumulation rates in the west, SR is even more restricted to low-elevation areas.
Two positive anomalies for both surface melt and surface runoff, widespread through all regions, took place in 1993 and 2006. They coincide with the highest melt index and melt extent found by Tedesco & Monaghan (Reference Tedesco and Monaghan2009) for the entire Antarctic continent and are associated with the Southern Annular Mode (SAM) and Southern Oscillation Index (SOI) negative anomalies. It demonstrates that the AP surface melt and surface runoff are closely linked to global climatic patterns and oscillations but remain poorly understood on the local spatial scales.
Conclusions
Using the global reanalysis ERAI, the present researchers ran a PDD-based model for surface melt and surface runoff estimation, calibrated and validated with local measurements. The results are in good agreement with previous local, regional and continental studies and can be used to reasonably investigate the combined influence of latitudinal, longitudinal and elevation differences along the Antarctic Peninsula on the spatiotemporal variability of surface melt and surface runoff. The ERAI-derived distributed PDDs presented here are strongly correlated with weather station-derived PDDs. By applying a melt factor of 5.4 mm PDD-1, good agreement with local measurements of SM was found.
The entire 1981−2014 period averages are 75Gt (surface melt) and 9Gt (surface runoff), with a very high interannual variability (σ M =54Gt and σ Q =8Gt). Although previous studies have found an exponential relation between PDDs and annual mean air temperature, the present work observed that the PDDs' interannual variability (and thereafter the amount of surface melt) is not well represented by the mean air temperature. Maximum values occurred in 1985 for surface melt (129Gt) and 1993 for surface runoff (40Gt). Since 2008, both variables have shown persistently negative anomalies, with minimum absolute values occurring in 2014 (surface melt=26Gt, surface runoff=0.37Gt). Nevertheless, no statistically significant temporal trends are present in the time series discussed here.
From 1996–2002, persistent negative anomalies are observed only for surface runoff. During this period, surface melt retention was always above 95% on the post-2002 floating ice areas. It is suggested that this persistent high retention may be linked to the Larsen B breakup mechanisms that took place in 2002. The post-2002 floating ice areas are responsible for a surface runoff average of 68% on the AP, highlighting their key importance to coastal hydrography as a freshwater source.
By dividing the AP into 11 regions, it can be observed that both the largest specific surface melt and runoff occur in the South Shetlands Islands, followed by the vicinity of the Larsen B Ice Shelf. This finding reveals the importance of the elevation and climatic differences between the western and eastern sides of the AP. Due to the Larsen C Ice Shelf, the east side has a flatter elevation profile, but generally also a smaller altitudinal decrease rate in temperature, and lower accumulation rates than the west side, resulting in higher specific surface melt and runoff in the east, despite the higher warming trends recorded in weather stations located in the west.
The large discrepancies among studies considering the same area and period indicate that further efforts are required to provide a better spatial distribution of field-measured surface mass balance data. The same recommendation can be made for weather stations. The discrepancies are possibly due to the high complexity of global, regional and local climatic forcing combined with the equally complex internal structure of the snow and ice layers.
Acknowledgements
Funding for this work was provided by EU FP7-PEOPLE-2012-IRSES IMCONet grant 318718, the Brazilian Coordination for the Improvement of Higher Education Personnel (CAPES) and the Brazilian National Institute of Science and Technology of the Cryosphere (INCT da Criosfera). M.B. was also funded via the German Research Foundation (grants BR2105/9-1 and BR2105/13-1). The RAMP DEM was kindly provided by NSIDC. Weather stations data were downloaded from NOAA at the website http://www.ncdc.noaa.gov/. ERA-Interim data were provided by the ECMWF. The surface mass balance data of Vega Island were kindly provided by Dr Sebastian Marinsek. The World Glacier Monitoring System distributes the surface mass balance data of Livingston Island, which are provided by Prof Dr Francisco Navarro. The authors would like to thank Prof Dr Francisco Navarro and the anonymous reviewer for their detailed revision and suggestions to improve the quality of this work.
Author contribution
J.C. conducted all data analyses and wrote major parts of the paper. M.B. and B.M. provided in situ reference data. N.E.B. provided the QuickSCAT-derived melt maps. J.C.S., A.B.S. and W.C.M. contributed to the data analysis and writing of the paper. J.C., J.A.-N. and M.B. jointly led the study. All authors contributed to and revised the manuscript. B.M. was funded by Russian Antarctic Expedition and IAEA project INT5153.
Details of data deposit
The new data presented in this paper will be archived at the PANGAEA data server (www.pangaea.de). PANGAEA is an information system hosted by the Alfred Wegener Institute, Helmholtz Center for Polar and Marine Research (AWI) and the Center for Marine Environmental Sciences, University of Bremen (MARUM).
