Deborah Verfaillie
email: deborah.verfaillie@ujf-grenoble.fr

Vincent Favier

Marie Dumont

Vincent Jomelli

Adrien Gilbert Gilbert

Daniel Brunstein

Hubert Gallée

Vincent Rinterknecht

Martin Menegoz

Yves Frenot

Recent glacier decline in the Kerguelen Islands (49°S, 69°E) derived from modeling, field observations, and satellite data

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Southern subpolar regions, located between 45 and 55 • S, play a key role in global atmospheric circulation [e.g., [START_REF] Garreaud | Present-day South American climate[END_REF][START_REF] Thompson | Signatures of the Antarctic ozone hole in Southern Hemisphere surface climate change[END_REF][START_REF] Purich | Autumn precipitation trends over Southern Hemisphere midlatitudes as simulated by CMIP5 models[END_REF] because energy from the tropics is mixed with polar cold air masses, leading to significant cyclonic activity in the region. Due to the association of low temperature with high humidity transported by the westerlies, large ice caps can form in Patagonia, New Zealand, and subantarctic islands despite the low elevation of their mountain ridges [e.g., [START_REF] Fitzharris | Glacier balance fluctuations and atmospheric circulation patterns over the Southern Alps, New Zealand[END_REF][START_REF] Takeuchi | Comparison of heat balance characteristics at five glaciers in the Southern Hemisphere, Global Planet[END_REF][START_REF] Fitzharris | Glacier balance fluctuations and atmospheric circulation patterns over the Southern Alps, New Zealand[END_REF][START_REF] Garreaud | Present-day South American climate[END_REF]. Similar to many other regions, glaciers in the southern subpolar regions are currently retreating. While glaciers in New Zealand retreated less and even advanced slightly from 1976 to the early 2000s, but have rapidly retreated since [START_REF] Fitzharris | Glacier balance fluctuations and atmospheric circulation patterns over the Southern Alps, New Zealand[END_REF], glaciers in Patagonia retreated strongly throughout the period, with a surface area loss between 0.4% and 36% since 1945 [START_REF] Lopez | A regional view of fluctuations in glacier length in southern South America, Global Planet[END_REF][START_REF] Melkonian | Satellite-derived volume loss rates and glacier speeds for the Cordillera Darwin Icefield, Chile[END_REF], and an acceleration of the retreat after the 1990s [START_REF] Rignot | Contribution of the Patagonia Icefields of South America to sea level rise[END_REF]Willis et al., 2012aWillis et al., , 2012b]].

The rapid response of glaciers to global warming has been widely studied by remote sensing, mainly on the Patagonian Ice Fields, where the detailed mass loss data were obtained from satellite imagery [e.g., [START_REF] Rignot | Contribution of the Patagonia Icefields of South America to sea level rise[END_REF]Willis et al., 2012aWillis et al., , 2012b]]. Glaciological modeling studies have also been performed in Patagonia [START_REF] Schaefer | Modeling past and future surface mass balance of the Northern Patagonia Icefield[END_REF][START_REF] Lenaerts | Extreme precipitation and climate gradients in Patagonia revealed by high-resolution regional atmospheric climate modelling[END_REF], but the models were forced by reanalysis data that contain biases (mainly concerning humidity) and did not include certain relevant parameters, such as calving processes. A limited amount of field data are available on the mass balance of southern subpolar glaciers. Indeed, long-term surface mass balance measurement programs have rarely continued for periods of more than a few years (with some exceptions such as Brewster glacier (New Zealand) since 2005 [START_REF] Anderson | Climate sensitivity of a high-precipitation glacier in New Zealand[END_REF]), and as a result, our knowledge of glaciological processes in the area is rather poor. In addition, the surface energy balance has only been analyzed for periods of less than 1 year [e.g., Hay and [START_REF] Hay | The synoptic climatology of ablation on a New Zealand glacier[END_REF][START_REF] Neale | Energy balance and synoptic climatology of a melting snowpack in the Southern Alps, New Zealand[END_REF][START_REF] Takeuchi | Comparison of heat balance characteristics at five glaciers in the Southern Hemisphere, Global Planet[END_REF][START_REF] Schneider | Energy balance in the ablation zone during the summer season at the Gran Campo Nevado Ice Cap in the Southern Andes[END_REF], and data are too scarce to clearly depict the main climatic changes that have occurred in the subpolar regions in recent decades.

The Kerguelen Islands (49 • S, 69 • E) are located in the southern Indian Ocean (Figure 1a) and form a unique subpolar observation site due to the scarcity of glacier-covered zones in this area. As these glaciers are located at low altitudes and on islands, they are directly influenced by oceanic and associated atmospheric variations [e.g., Poggi, 1977aPoggi, , 1977b;;[START_REF] Vallon | Glaciologie à Kerguelen[END_REF], and they have been subject to marked fluctuations in the last 50 years, while their state had been relatively stable between the 1800s and the 1960s [START_REF] Frenot | Azorella selago Hook. Used to estimate glacier fluctuations and climatic history in the Kerguelen Islands over the last two centuries[END_REF]. After the 1960s, the Cook Ice Cap retreated extremely rapidly, losing 20% of its surface area in 40 years [START_REF] Berthier | Ice wastage on the Kerguelen Islands (49 • S, 69 • E) between 1963 and 2006[END_REF]. This acceleration has generally been attributed to warming, while fluctuating precipitation has generally been assumed to be of less importance [e.g., [START_REF] Frenot | Azorella selago Hook. Used to estimate glacier fluctuations and climatic history in the Kerguelen Islands over the last two centuries[END_REF][START_REF] Frenot | Datation of some Holocene peat sediments and glacier fluctuations in the Kerguelen Islands[END_REF][START_REF] Berthier | Ice wastage on the Kerguelen Islands (49 • S, 69 • E) between 1963 and 2006[END_REF]. However, due to our lack of knowledge of the relationship between the local climate, the glaciers, and ice dynamics in the area, until recently it was difficult to conclude whether the recent drying trend had a serious impact on glacier wastage.

As part of a 4 year glaciological program in the Kerguelen Islands, in 2010, a meteorological, glaciological, and hydrological network was set up on and around Ampere Glacier-an outlet glacier of the Cook Ice Cap (Figure 1a)-to detect, monitor, and understand climate and mass balance variability and change in the glacial environment of the ice cap. We also used data from previous short-term mass balance studies [Vallon, 1977a[Vallon, , 1977b[START_REF] Vallon | Glaciologie à Kerguelen[END_REF]] and short-term energy balance studies [e.g., Poggi, 1977aPoggi, , 1977b] ] conducted on Ampere Glacier between 1970 and 1974 for comparison with our new results.

The main objective of the present study was to document the recent ice cap wastage of this poorly studied subpolar area. Additionally, a first comparison with other subantarctic glacial sites such as Patagonia and New Zealand was also carried out. Ultimately, this improved our knowledge of the relationship between the glaciers and local (and regional) climate.

The field data (section 2.1) analyzed in this study are the first step in understanding the recent rapid retreat of the cryosphere on Kerguelen and were used to validate the surface mass balance and surface energy Journal of Geophysical Research: Earth Surface 10.1002/2014JF003329 balance models. But in addition, collecting spatially distributed data on the surface characteristics and ablation from satellite images (section 2.2) was crucial to validate our models.

Following [START_REF] Lopez | Snow cover monitoring in the Northern Patagonia Icefield using MODIS satellite images (2000-2006[END_REF], we produced our own images of the ice and snow cover from MODIS Level 1B images. We applied the method developed by [START_REF] Sirguey | Subpixel monitoring of the seasonal snow cover with MODIS at 250 m spatial resolution in the Southern Alps of New Zealand: Methodology and accuracy assessment[END_REF] and [START_REF] Dumont | Linking glacier annual mass balance and glacier albedo retrieved from MODIS data[END_REF] to produce surface albedo maps from MODIS Level 1B images of the Cook Ice Cap to estimate the transient snow line and assess changes in the extent of the accumulation area of the ice cap since 2000. In addition, MODIS albedo values measured at the location of stakes enabled computation of the local surface mass balance through energy balance estimates from a surface energy balance model, as well as comparison with surface mass balance values measured in the year 2011. To identify the causes and processes involved in the retreat of the ice cap, the third method we used was analyzing climatic data (section 2.3) over Kerguelen derived from measurements and reanalyses since the 1950s. Finally, we obtained an additional estimate of surface mass balance at the ice cap scale by implementing a simple positive degree-day model based on temperature and precipitation measured at the scientific station of Port-aux-Français. The model was validated using surface energy and mass balance measurements. The modeling results were then compared with previous estimates derived from field observations and satellite imagery.

Data

Field Data on and Around the Glacier

Field data from a meteorological and glaciological network set up on and around Ampere Glacier (Figure 1b) were used in this study.

These data were obtained from field campaigns conducted in 2010-2011 and 2011-2012, during which ablation stakes were measured to assess the glacier surface mass balance (SMB). SMB is the sum of the accumulation at the surface of the ice sheet (through precipitation, hoar formation, and wind deposition) and surface ablation (through sublimation, melting, and wind scouring). Twenty ablation stakes were set up in December 2010 and January 2011 and measured 1 year later in December 2011 and January 2012 (stake height and location were measured). The stakes were separated into two categories (Figure 1b): those located in the regular ablation area (zone 1) and those located at specific points where snow overaccumulates due to wind deposition (zone 2). An automatic weather station (AWS) was also set up on a rock outcrop close to the glacier (Figure 1b, La Mortadelle, 200 m above sea level (asl)) to record local meteorological conditions. The characteristics of the different AWS sensors are listed in Table 1. Thirty minute averaged data on air temperature and humidity, wind speed and direction, precipitation, and incoming and outgoing shortwave and longwave radiation (SW and LW, respectively) were measured at the AWS. Additional surface short-term heat budget measurements were recorded by a second AWS (Figure 1b, Ampere, 280 m asl) placed on the glacier surface during the field campaigns. A stake farm was also installed around this AWS. It consisted of four stakes that were measured daily and used for calibration of the SMB model (section 3.5).

Additionally, data from previous surface energy balance (SEB) studies [Poggi, 1977a] and SMB studies using stakes [Vallon, 1977a[Vallon, , 1977b[START_REF] Vallon | Glaciologie à Kerguelen[END_REF], conducted on Ampere Glacier in the 1970s (the positions of the stakes are shown in Figure 1b) were also analyzed for comparison.

Journal of Geophysical Research: Earth Surface 10.1002/2014JF003329

Satellite Images 2.2.1. MODIS Sensor

The MODIS sensor is a multispectral radiometer on board the EOS/AM-1 Terra platform, working in 36 bands from the visible to the thermal infrared. It produces images of the Earth's surface since February 2000, at three spatial resolutions: 250 m (bands 1 and 2), 500 m (bands 3 to 7), and 1 km (all remaining bands). We used MODIS images for our study because of its high temporal resolution (global coverage in less than 2 days and more than one image per day over Kerguelen), which is required to study changes in the seasonal snow line on the Cook Ice Cap. The different MODIS products used in this study are described in the following sections.

MODIS Level 1B Products

For the production of albedo maps, we used MODIS Level 1B Swath images (MOD02QKM, MOD02HKM, MOD021KM, and MOD03) distributed by the Level 1 and Atmosphere Archive and Distribution System (LAADS, http://ladsweb.nascom.nasa.gov/). We preferred these images to MOD10 and MCD43 products because of ( 1) their finer resolution (250 m for MOD02QKM);

(2) the possibility to retrieve different types of albedo from them (as explained in section 3.1); (3) their high accuracy (root-mean-square deviation of about 0.05 between MODIS-retrieved albedo and field-measured albedo as shown in [START_REF] Dumont | Linking glacier annual mass balance and glacier albedo retrieved from MODIS data[END_REF]); and (4) a misinterpretation of Kerguelen snow surfaces in MOD10, which, in those images, are mapped as clouds (also the case in some regions of New Zealand [START_REF] Hall | MODIS snow-cover products[END_REF][START_REF] Sirguey | Subpixel monitoring of the seasonal snow cover with MODIS at 250 m spatial resolution in the Southern Alps of New Zealand: Methodology and accuracy assessment[END_REF] and in Patagonia [START_REF] Lopez | Snow cover monitoring in the Northern Patagonia Icefield using MODIS satellite images (2000-2006[END_REF]). MODIS Level 1B Swath products contain radiance measurements at the top of atmosphere in different bands and a geolocation file providing geometry and topography parameters (such as solar and sensor azimuths and zeniths).

The Cook Ice Cap behaves similarly to alpine glaciers, with an ablation season mainly during (austral) summer and an accumulation season mainly during (austral) winter [Vallon, 1977a]. Consequently, we analyzed images in three different periods between 2000 and 2012: (1) November-December, which corresponds to the beginning of the ablation season;

(2) March-April, which is the end of the ablation season; and (3) July-August, in the middle of winter, hence the accumulation season.

Those periods were chosen to evaluate the seasonal variations in the location of the snow line and its highest position (in March-April), as well as any differences between 2000 and the present. A complementary analysis of the variations in the snow line of the Ampere Glacier was also carried out. We only used images of the Cook Ice Cap and the Ampere Glacier that were completely clear, and thus suitable for glacier monitoring, which came to a total of only 57 images out of the 8700 images available. In addition, between December 2010 and December 2011, one to three images per month of the Ampere Glacier (but not always the whole of Cook Ice Cap) were used to estimate variations in the surface albedo and hence to calculate the SW radiation budget for 2011 (see section 3.4).

Fine-Resolution Images

Additional Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and Landsat images were used to draw a mask of the ice cap for each albedo map (using the image taken on the date closest to the date of the albedo map). To validate the albedo retrieval method (section 3.1) and the threshold, we chose for the albedo of snow (section 3.2), and we compared one of our albedo maps with an ASTER image taken on the exact same date. ASTER and Landsat images have a much finer spatial resolution (15 m and 30 m, respectively) than MODIS images but are only rarely available. Complementary fine-resolution images used in this study are described in Table 2. Two digital elevation models (DEM) were also used to assess the elevation of the ice cap in 2000 and 2009: one generated from the February 2000 Shuttle Radar Topographic Mission (SRTM) [START_REF] Farr | The Shuttle Radar Topography Mission[END_REF] and the other from Satellite Pour l'Observation de la Terre (SPOT-5) optical stereo imagery acquired in December 2009 [START_REF] Korona | SPIRIT. SPOT 5 stereoscopic survey of Polar Ice: Reference images and topographies during the fourth international polar year (2007-2009[END_REF]. Surface elevation was used to determine snow line elevations since 2000 (by considering the DEM with the date closest to the acquisition of each MODIS image) and to estimate ice cap hypsometry in 2009.

Meteorological and Climatic Data

Direct meteorological measurements have been collected routinely by Météo France at the station of Port-aux-Français (29 m asl; see Figure 1a) since 1950. These data were used to get a clear view of local climatic variations and for long-term SMB modeling. Measurements were also compared to model results using complementary reanalysis data to see whether variations were regional or only local.

Journal of Geophysical Research: Earth Surface 10.1002/2014JF003329 We used ERA-40 and ERA-Interim data [START_REF] Uppala | The ERA-40 re-analysis[END_REF][START_REF] Dee | The ERA-Interim reanalysis: Configuration and performance of the data assimilation system[END_REF] from the European Centre for Medium-Range Weather Forecasts (ECMWF) and Reanalysis 1 and 2 data [START_REF] Kalnay | The NCEP/NCAR 40-year reanalysis project[END_REF][START_REF] Kanamitsu | NCEP-DOE AMIP-II Reanalysis (R-2)[END_REF] from the National Centers for Environmental Prediction, National Center for Atmospheric Research, and Department of Energy (NCEP-NCAR, NCEP-DOE). Reanalyses correspond to outputs from models using data assimilation, i.e., combination-during a given cycle-of observations (from radiosondes, buoys, satellites, etc.) and prior forecast model information, to obtain a physically coherent estimate of the state of the atmosphere that is as close to reality as possible. Each reanalysis is produced with a single version of a data assimilation system. Surface and radiosonde observations performed at Port-aux-Français have been assimilated in the reanalyses since 1967 and 1968 respectively for ERA and 1967 and 1970 respectively for NCEP.

As grids and resolutions differ between models (i.e., ∼210 km for NCEP 1 and NCEP 2, ∼125 km for ERA-40, and ∼80 km for ERA-Interim), reanalysis data for Kerguelen were retrieved for different grid cell locations depending on the reanalysis concerned. Grid cells used in this study are centered on 50 • S, 70

Methods

Albedo Retrieval Method

The MODIS sensor provides values of radiance measured at the top of the atmosphere. To obtain surface albedo, several algorithms were developed by [START_REF] Sirguey | Subpixel monitoring of the seasonal snow cover with MODIS at 250 m spatial resolution in the Southern Alps of New Zealand: Methodology and accuracy assessment[END_REF] and [START_REF] Dumont | Linking glacier annual mass balance and glacier albedo retrieved from MODIS data[END_REF]. A first algorithm was performed to retrieve the subpixel snow cover and hemispherical-conical reflectance for seven spectral bands from MODIS images. This algorithm accounts for the correction of atmospheric effects and the multiple reflections that occur in rugged terrain. The pixels identified as snow or ice are then corrected using another algorithm to obtain directional-hemispherical reflectance. A third and last algorithm transforms directional-hemispherical reflectance into bihemispherical reflectance and then bolometric albedo. Blue sky albedo is obtained as well as white sky albedo, the latter being the albedo of the same surface under pure diffuse illumination. The albedo we used in this study was the white sky albedo, as it makes it possible to compare albedo maps obtained at different periods of the year when the solar zenith angle is not the same. Albedo maps (250 m resolution) were finally obtained from each processed MODIS image.

Determination of the Transient Snow Line and Accumulation Area

The equilibrium line altitude (ELA) can generally be approximated to the altitude of the snow line (separation between the area covered by snow and the area where ice is visible) at the end of the ablation season [e.g., [START_REF] Rabatel | Can the snowline be used as an indicator of the equilibrium line and mass balance for glaciers in the outer tropics?[END_REF]. Figure 2 shows the comparison between an albedo map derived from MODIS and an ASTER image taken on the same day, which helped us determine the albedo threshold between ice and snow. Figure 2b shows that an albedo threshold of 0.65 appears to mark the limit between ice and recent bright snow (area above the continuous black line), thus marking the position of the transient snow line (i.e., the line separating snow surfaces from ice or firn surfaces at any instant [START_REF] Cogley | Glossary of Glacier Mass Balance and Related Terms[END_REF]). However, especially in the NW part of the ice cap and on the ablation zone of Ampere Glacier (Figure 2b), darker areas appear between an albedo threshold of 0.55 (dotted black line) and 0.65 (continuous black line). Upon closer look (Figure 2c), these areas seem to correspond to fragmented fresh snow deposits where snow accumulated in topographic lows or at the lee of crests. To account for the presence of these fragmented snow deposits, we chose a snow line albedo threshold of 0.6 ± 0.05.

This threshold is slightly higher than the commonly measured snow line albedo value. For example, the maximum albedo value of the ablation zone on Saint-Sorlin glacier (French Alps), derived from AWS measurements, terrestrial photographs, and MODIS images is 0.5 [START_REF] Dumont | Monitoring spatial and temporal variations of surface albedo on Saint Sorlin Glacier (French Alps) using terrestrial photography[END_REF][START_REF] Dumont | Linking glacier annual mass balance and glacier albedo retrieved from MODIS data[END_REF]. There are three possible reasons for this difference in albedo: (1) fewer impurities in Kerguelen air masses than in the French Alps;

(2) the relative small amount of rock outcrops (nunataks) on Ampere Glacier compared to alpine glaciers, which are generally surrounded by rock faces; and (3) very frequent snowfalls all year round on Kerguelen, which cause frequent increases in the snow albedo. Given the difference in resolution between the MODIS and ASTER images, the good visual correspondence between snow areas in the two images gives us confidence in our method of albedo retrieval and our choice of albedo threshold.

Pixels with an albedo exceeding 0.6 ± 0.05 located inside the ice cap mask were counted to assess the snow-covered area for each albedo map. The total area of the ice cap was also estimated for each albedo map (albedo > 0.3). The ratio of snow-covered areas to areas not covered with snow in the Cook Ice Cap was then calculated for each image.

Surface Energy Balance Model

Point SEB computations were performed using data from La Mortadelle and Ampere AWS and compared to estimates of point energy budget in section 3. surface of the glacier (Figure 1b) to measure albedo, air temperature and humidity, and wind speed. The energy stored in the top layers of the glacier was calculated as (in W m -2 , fluxes toward the surface are positive) [START_REF] Oke | Boundary Layer Climates[END_REF][START_REF] Favier | Modeling the mass and surface heat budgets in a coastal blue ice area of Adelie Land, Antarctica[END_REF]]

SW i + SW r + LW i -((1 -𝜖)LW i + 𝜖𝜎T 4 s ) + LE + H = F surface , ( 1
)
where SW i is incoming shortwave radiation and SW r is reflected shortwave radiation that were directly measured on the glacier; LW i is incoming longwave radiation that was measured at La Mortadelle AWS, and the term in brackets is outgoing longwave radiation LW r that was deduced form the modeled surface temperature T s and Stefan-Boltzmann's equation. Surface emissivity is 𝜖 = 0.99 and the Stefan-Boltzmann constant is 𝜎 = 5.67 10 -8 W m 2 K -4 . H and LE are turbulent and sensible heat fluxes, respectively, which were computed using the bulk method including stability corrections based on the bulk Richardson number. Surface roughness lengths were taken from the literature [START_REF] Favier | Modeling the mass and surface heat budgets in a coastal blue ice area of Adelie Land, Antarctica[END_REF], to get agreement between the ablation measured at the stakes and modeled ablation. The heat flux supplied by precipitation was not taken into account because during precipitation, the temperature is always close to 0 • C.

F surface is the energy available at the surface. Some shortwave radiation is not available for warming/cooling processes at the surface, or for melting, because the shortwave flux partly penetrates the ice. Hence, F surface is separated into two terms:

F surface = G 0 + (1 -a)SW, ( 2
)
where G 0 is the energy excess or deficit at the surface. In this equation, a is the fractional amount of shortwave radiation that is absorbed in the top layer of the model (at the surface).

When the surface temperature is 0 • C, positive G 0 values represent the energy available for melt. Otherwise, this amount is used to cool/warm the frozen surface and underlying snow/ice. For this task, heat conduction within the snow or ice was computed by solving the thermal diffusion equation according to an explicit scheme to a depth of 2 m, with a 5 cm grid resolution and a 20 s time step.

Running a full surface energy balance model was not possible over the last decades, because (1) accurate measurements of the energy fluxes at the glacier surface were only available from 21 December 2010 to 4 January 2011 and from 14 December 2011 to 30 December 2011, (2) long-term radiation (only shortwave) and wind speed measurements are available at Kerguelen since the 1990s only, and (3) the reanalyzed (ERA-40 and NCEP 1) radiation and cloudiness data at Kerguelen do not significantly correlate with measurements performed at Kerguelen since 1990 (not shown). As Ampere AWS data were only available during field campaigns, annual outputs from this SEB model have not been validated yet, but we are confident in the estimates provided by the model (presented below). This confidence is reinforced by the good correspondence between surface height estimates using the SEB approach and using the PDD approach (section 3.5).

Estimation of Point SMB Using MODIS Albedo and SEB Estimates

Based on MODIS albedo estimates, we calculated point SMB at the location of the stakes on the basis of estimated solid precipitation and a simple computation of ablation, in which the SW radiation budget is calculated using MODIS albedo values (section 3.2). In addition, we estimated SMB for a pixel located in the accumulation zone (see Figure 1b), at 990 m asl, for which all the images used for computation of SMB at the ablation stakes were clear. Our method is similar to the one used by [START_REF] Greuell | Assessment of interannual variations in the surface mass balance of 18 Svalbard glaciers from the Moderate Resolution Imaging Spectroradiometer/Terra albedo product[END_REF], except that we calculated point SMB values instead of glacier-averaged values.

The annual accumulation term at the location of each stake can be estimated as

Accumulation = n ∑ i=1 P S , ( 3
)
where P S is the amount of solid precipitation (in mm), occurring during events i = 1 to n. Hoar formation and wind deposition are not taken into account. [START_REF] Jóhannesson | Degree-day glacier mass-balance modelling with applications to glaciers in Iceland, Norway and Greenland[END_REF] and is equal to 0 otherwise. T air is the temperature from the Mortadelle AWS data; z is stake altitude (in m asl); 𝛾 is the lapse rate, used to account for the effect of altitude on precipitation phase; and P is the 30 min precipitation amount on Ampere Glacier estimated from measured precipitation at Port-aux-Français extrapolated

P S = P if T air + 𝛾 × z ≤ 1 • C [e.g.,
Journal of Geophysical Research: Earth Surface 10.1002/2014JF003329

to the glacier using precipitation distribution given by the regional model MAR (Modèle Atmosphérique Régional) [START_REF] Gallée | Development of a three-dimensional meso-𝛾 primitive equation model: Katabatic winds simulation in the area of Terra Nova Bay, Antarctica[END_REF].

The MAR model is a hydrostatic climate model using the atmospheric scheme described in [START_REF] Gallée | Development of a three-dimensional meso-𝛾 primitive equation model: Katabatic winds simulation in the area of Terra Nova Bay, Antarctica[END_REF], coupled to a physically based model of the snow pack [START_REF] Gallée | Air-snow interactions and the surface energy and mass balance over the melting zone of west Greenland during the Greenland Ice Margin Experiment[END_REF][START_REF] Gallée | Development of a three-dimensional meso-𝛾 primitive equation model: Katabatic winds simulation in the area of Terra Nova Bay, Antarctica[END_REF]. The simulation was run at 10 km resolution on a stereographic grid for the year 2011. The model was forced with ERA-Interim reanalysis data (see section 2.3). Results from the MAR model give a lapse rate of -8.0 • C/km along the azimuth of Ampere Glacier (southeast) and of -8.7 • C/km eastward.

For Ampere Glacier, we verified that the computed lapse rate value agreed with the lapse rate computed with meteorological stations at La Mortadelle and on Ampere Glacier (-9.1 • C/km in summer). For glaciers oriented to the east, the lapse rate was validated with a comparison between La Mortadelle and Port-aux-Français, giving a lapse rate value of -8.6 • C/km in 2011 and 2012. For the calculation of P, briefly, the relationship between modeled elevation and modeled total precipitation in the MAR model is

P(z) = P PAF (2.79 z + 783)∕783, ( 4
)
where P PAF is the 30 min precipitation amount measured at Port-aux-Français. This relationship was used to calculate total precipitation and then solid precipitation at stake altitude on Ampere Glacier.

Because of the high amount of humidity in the air, sublimation is low on the Cook Ice Cap [e.g., Poggi, 1977a]. Since melting amounts are large (7 to 10 m water equivalent (w.e.) per year in the ablation area), the contribution of sublimation to the ablation term can be ignored (see section 3.3). Similarly, we ignored the contribution of wind scouring. In addition, conduction through ice and snow (G) can also be ignored because the ice cap is temperate. Consequently, the annual ablation term can be estimated as

Ablation ≈ Melting = (SW + LW + LE + H + Pr) 1 L f , ( 5
)
where SW is the shortwave radiation budget, LW is the longwave radiation budget, LE and H represent the turbulent latent and sensible heat contributions, and Pr is the contribution of precipitation. Pr is small compared to other terms [Poggi, 1977a] and can thus be ignored. Finally, L f is the latent heat of fusion (= 0.334 × 10 6 J kg -1 ). Contributions are annual energy amounts expressed in W m -2 a -1 .

We calculated the SW contribution to SMB at each stake location on Ampere Glacier as

SW ∼ 1 n n ∑ m=1 [(1 -ᾱm )SW m i ], (6)
where n = 12 (the number of months in a year), ᾱm is the mean albedo obtained during month m from available clear-sky images (one to three images per month), and SW m i is the monthly accumulated incoming SW radiation measured at La Mortadelle AWS during month m. For months with no clear-sky images, we used the albedo value from the closest date.

The incoming longwave radiation (LW i ) on the glacier was estimated from La Mortadelle AWS measurements assuming that there are no major variation in the mean flux (mean annual flux is 320 W m -2 ) between 200 m asl and 600 m asl. The outgoing longwave radiation flux was estimated using the Stefan-Boltzmann law, yielding a value of -311 ± 5 W m -2 assuming a mean surface temperature between -2 • C and 0 • C. This surface temperature range was chosen given that (1) the Mortadelle AWS monthly air temperature in winter is usually above 0 • C and (2) some preliminary point SEB simulations (section 3.3) yielded a mean annual surface temperature of -1 • C. These assumptions suggest a mean annual longwave radiation budget of 9 ± 5 W m -2 . LE and H were not directly measured on our AWS. As the state and temperature of the glacier surface are not expected to vary much between 200 m asl and 600 m asl because of almost continuous melt rates, a fixed value for turbulent heat fluxes was chosen that does not depend on elevation. A rough estimate of the turbulent heat contribution was first assessed from Poggi [1977a] (50 W m -2 ) and then compared to estimates from a point SEB model (34 W m -2 , section 3.3).

Uncertainties on SMB due to albedo were roughly estimated as the highest error introduced when we used different pixels adjacent to the one in which the stake was located. This yielded uncertainties of 6% for the accumulation point, ±10% for stakes located in zone 1 and between ±10% and ±31% for stakes located in zone 2. In addition, because MODIS albedo estimates were only considered under clear-sky conditions, Thin black and gray lines with squares are the ice and snow surface height (respectively) computed with the PDD approach. Black and gray circles are the mean measured ice and snow surface height (respectively) at a stake farm (four stakes) located around the Ampere AWS on the ablation zone (280 m asl; see Figure 1b). The red curve is the 30 min surface elevation measured with an acoustic gauge at the same location. Error bars are ±20 mm, which represent the mean standard deviation between stake measurements. a slight underestimation of the mean annual albedo occurs. Indeed, (1) albedo in clear-sky conditions is lower than when some cloud cover is present because there is less absorption of near-infrared radiations (i.e., a part of the solar spectrum for which the albedo is lower) and (2) fresh snow may already have melted in clear-sky images. A rough estimate of the resulting uncertainty suggests that albedo values should be increased by +0.05, leading to a systematic bias on SMB values of +0.5 m w.e. a -1 .

Positive Degree-Day Model

We computed the snow line and SMB with a simple positive degree-day (PDD) model based on [START_REF] Hock | A distributed temperature-index ice-and snowmelt model including potential direct solar radiation[END_REF] and fully described in [START_REF] Maisincho | On the interest of positive degree day models for mass balance modeling in the inner tropics[END_REF]. Solid precipitation on the ice cap was estimated using the same method as described in section 3.4, with precipitation data from Port-aux-Français available since 1951.

Figure 3 shows the comparison of computed and measured surface height on the ablation area (280 m asl) using meteorological data from Ampere AWS (and precipitation from La Mortadelle AWS), which was used to calibrate the ice and snow melting coefficients. Comparison with outputs from the SEB model is also shown, indicating a good agreement between surface height computed with the PDD model and with the SEB model. In December 2010 (Figure 3, top), a lot of minor snowfall events influenced measurements as well as modeling. As a result, the melting coefficient for ice was calibrated using measurements made in December 2011 (Figure 3, bottom). However, the lack of situations with long-term snow cover led us to use a snow coefficient from the literature, which is the average of values from [START_REF] Radic | Regionally differentiated contribution of mountain glaciers and ice caps to future sea-level rise[END_REF] (4.9 mm w.e. • C -1 d -1 ). For ice, the value in this study (7.4 mm w.e. • C -1 d -1 ) is close to the average value from [START_REF] Radic | Regionally differentiated contribution of mountain glaciers and ice caps to future sea-level rise[END_REF] (7.2 mm w.e. • C -1 d -1 ). Both values used for Ampere Glacier are also close to values determined empirically in [START_REF] Anderson | Past and future mass balance of 'Ka Roimata o Hine Hukatere' Franz Josef Glacier[END_REF] for Franz Josef Glacier (4.55 ± 0.50 mm w.e. • C -1 d -1 for snow and 7.17 ± 0.20 mm w.e. • C -1 d -1 for ice).

Modeling should reproduce the reduction in size of the area covered by firn (not visible in Figure 2 but noticed during field campaigns). To account for this process, an initial firn thickness from 1951 was used on which a PDD run was performed over 20 years (spin-up) to assess the initial firn cover in 1971, and the PDD run was then extended to the present. We assumed different cases with firn layers ranging from 5 to 15 m w.e., which are discussed in section 4.1.

In addition, a second run was performed without an initial firn cover, assuming that the firn layer at the beginning of the 1970s resulted only from the mean annual accumulation and ablation from the previous years.

Results

PDD Validation Using Field Data

SMB at the ice cap scale was modeled with a simple PDD model based on meteorological measurements (section 3.5). Results are shown in Figure 4 for comparison with stake measurements (presented in section 4.3). Our PDD model results are in agreement with SMB measurements (Figure 4) for the recent period (medium correlation (R) of 0.58 discussed below, p value <0.001 according to a Student's t test) and the 1970s (excellent correlation of 0.99, p value <0.001). This shows that the PDD model initiated with a 1b) plotted against surface elevation (m above sea level). Values obtained between 1970 and 1974 by Vallon [1977aVallon [ , 1977b[START_REF] Vallon | Glaciologie à Kerguelen[END_REF] are also given for the purpose of comparison. Lines represent ice cap SMB derived from our PDD model, for the early 1970s (red line, average of 1970-1974) and 2011 (blue lines). The PDD SMB for 2011 is calculated from January to January (light blue line) and from May to May (dark blue line). The dotted lines represent the version of the model with no initial firn cover, while the solid lines represent the version with an initial firn cover. Uncertainty intervals (±10 m w.e.) in the initial firn accumulation are also given. 20 year spin-up and a firn cover is capable of reproducing the SMB distribution with elevation and the SMB changes between the 1970s and 2011. However, Figure 4 shows that the initial amount of firn does not have much impact on modeled SMB. Indeed, a slight difference in SMB (only concerning elevations close to the firn and snow lines) can be observed between 5 m and 15 m firn experiments, but higher amounts of initial firn yielded the same final SMB as the 15 m firn experiment.

The PDD model shows that the snow line was above 750 m asl in 2011 (this was confirmed by field observations in 2012 and 2013), whereas it was located around 650 m asl in the early 1970s. Even if the snow line elevation in 2011 was determined only using data from 1 year, the PDD model shows that it was similar between 2008 and 2012, suggesting that this difference between 1970s and 2011 makes sense.

Recent Variations in the Snow Line and Accumulation Area Derived From MODIS Albedo Maps

Figures 5b and5c show changes in the ratio between areas covered by snow (albedo > 0.6 ± 0.05) on the ice cap and the entire ice cap surface between 2001 and 2012 calculated from albedo maps (section 3.2). This ratio may be considered equivalent to the accumulation area ratio (AAR) at the end of the ablation season, if we assume that the snow line elevation is similar to the ELA at the end of summer. In Figure 5, we use the term "AAR, " even though we acknowledge that winter values of this ratio are different from the strict sense AAR which should be based on annual values. Despite the marked scattering of the results and the limited number of good quality images, the ratio decreased between 2001 and 2012. A low snow line appeared on 5 April 2012 followed by a rapid increase that was revealed on 26 April 2012, suggesting that estimates may not reflect the highest annual snow line position, and the trends obtained may not be very accurate due to the small number of images available. However, the same trend was observed in each season, showing that the ratio had decreased by ∼20% at the end of the ablation season from 2001 to 2012 (which corresponds to a decrease of ∼90 km 2 in the snow-covered areas). Values for March-April 2008 and 2009 are, however, the highest of the March-April period.

Figure 6 shows two albedo maps obtained for the period March-April (end of the ablation season) 2001 and 2012. The difference in the snow line (dotted and continuous black lines) between the two images is striking, with an increase of about 100 m between those two periods. The images were not acquired at the exact same date but about 6 weeks apart due to the scarcity of clear-sky images. Thus, the disappearance of the snow cover could be partly linked to the time lag between the images, but PDD modeling showed that this reduction is consistent (section 4.4) and that the estimated snow line in 2001 was close to its maximum. Moreover, the MODIS images also showed a retreat of the ice cap between 2001 and 2012. Indeed, the disappearance of ice, unlike that of snow, does not depend on the time of year, and Figure 6 shows that some glacier snouts and termini have disappeared in 2012, especially on the eastern flank of the ice cap (sites are indicated by pink arrows). The glacier retreat was confirmed by comparing Landsat images acquired in 2001 and 2011 (not shown here). The point named accu corresponds to the x = y SMB value calculated with the albedo of a pixel situated in the accumulation zone of the glacier (990 m asl). Error bars were estimated by calculating the error introduced when selecting different pixels adjacent to the one in which the stake is located, plus the systematic bias introduced by using clear-sky images.

Current SMB and Snow Line Estimates From Stake Measurements and Comparison With the 1970s

extrapolation of stake measurements suggests that in 2011, the snow line was located at about 800-1000 m asl, in agreement with the transient snow line value obtained for Ampere Glacier from MODIS images in 2012 (Figure 6). However, the data display marked scattering caused by five values from stakes located at the border of the glacier where snow can accumulate (Figure 1b, zone 2).

Despite the two short observation periods, we were able to compare SMB values on Ampere Glacier between 1970 and 1974 [Vallon, 1977a[Vallon, , 1977b[START_REF] Vallon | Glaciologie à Kerguelen[END_REF] with our values for 2011 (Figure 4). The comparison revealed that the spatial distribution of SMB has changed significantly in 40 years: (1) SMB is more negative at the same elevation and (2) SMB gradient above 400 m asl is lower than in the early 1970s.

Comparison Between the PDD Snow Line and the MODIS Transient Snow Line

Figure 5a compares the modeled snow line and the transient snow line elevation obtained from MODIS images on Ampere Glacier. Modeled snow line elevations are underestimated compared to MODIS transient snow line elevations: the mean difference between the PDD-modeled snow line and MODIS-derived transient snow line (x PDD -x MODIS ) is -112 m for a transient snow line threshold of 0.6. However, despite its high variability, the modeled snow line reproduces well variations in the transient snow line from MODIS images (correlation of 0.72 for a transient snow line at an albedo of 0.6, p value <0.001). The annual maximum PDD-modeled snow line elevation (corresponding to the ELA) of the Ampere Glacier increased by ∼30 m between 2000 and 2012.

Information Provided by SEB Modeling and Albedo Estimates From MODIS

We calculated point SMB values at stake locations (see Figure 1b for locations) for 2011 and compared them to values measured at the stakes in 2011, as illustrated in Figure 7.

The poor quality of the information concerning turbulent heat fluxes and precipitation constitutes the main uncertainty in SEB and SMB computation. High LE + H values tend to decrease the local mass balance.

Conversely, high precipitation tends to limit negative SMB values. We observed that field SMB data cannot be reproduced if we assume a value of 50 W m -2 for turbulent heat fluxes (as measured by Poggi [1977a]) even if we assume a P twice as large as the one calculated from the relationship described in section 3.4 (Figure 7a). On the other hand, we observed that SMB is correctly computed assuming calculated precipitation and turbulent heat fluxes ∼30 W m -2 (Figure 7b). Preliminary results from our SEB model during 2011 yielded a similar mean value of point LE + H of 34 W m -2 . Assuming doubled precipitation to account for rain gauges catch deficiency in strong wind conditions [START_REF] Larson | Accuracy of precipitation measurement for hydrological modelling[END_REF] has little effect on final SMB estimates, except for high-altitude points (i.e., the "accu" point).

However, even for the case with turbulent heat fluxes of 30 W m -2 , points located in overaccumulation areas (zone 2) yield SMB estimates too negative compared to measurements (they fall under the x = y curve), while the ones located in zone 1 agree well with measurements. About 3 times the amount of P (which is mostly liquid) have to be considered for the computed SMB in zone 2 to agree with measurements, suggesting that they are located in large corniches. Finally, we observed that the point chosen at 990 m asl, where accumulation is expected to occur, displays a negative estimated SMB value (approximately -3200 mm w.e. for the case with turbulent heat fluxes of 30 W m -2 at an elevation of 990 m asl). In the case of LE + H = 30 W m -2 , accumulation calculated from the relationship described in section 3.4 (approximately +3350 mm w.e.) compensates the SW contribution to SMB (approximately -3140 mm w.e.). This demonstrates that LE + H and LW must decrease with elevation in order to have a slight accumulation at the top of the ice cap and that constant values used in the ablation zone are not appropriate at such high elevations.

Long-Term Variations in the Snow Line and Link With Climate Variability

Measured and reanalyzed temperature and precipitation data over Kerguelen since 1948 are presented in Figure 8.

Météo France measurements revealed a marked decrease in temperature of about -0.06 • C a -1 from 1952 to 1964, followed by sharp increase of about +0.06 • C a from 1964 to 1982. and measurements after 1982 revealed a slight increase (not significant at the 95% confidence level) in temperature of about +0.002 • C a -1 to +0.01 • C a -1 .

Precipitation measured at Port-aux-Français before 1973 presented a very different variability than precipitation measured after 1973. Measurements after 1973 displayed a significant (at the 95% confidence level) decrease in precipitation of about 25% in the last 40 years, slightly higher than in reanalysis data, which suggest a decrease of 5% to 20%. Analysis of daily precipitation revealed that this decrease is not linked to a decrease in the intensity of precipitation but rather to a decrease in the frequency of precipitation events (Table 3). However, it should be kept in mind that precipitation measurements using rain gauges are negatively biased in strong wind conditions [START_REF] Larson | Accuracy of precipitation measurement for hydrological modelling[END_REF], which is often the case on the Kerguelen Islands, implying that actual precipitation amounts are probably larger than measured but measurements were likely uniformly impacted during the last 60 years. [START_REF] Frenot | Azorella selago Hook. Used to estimate glacier fluctuations and climatic history in the Kerguelen Islands over the last two centuries[END_REF]). Considering unchanged precipitation amounts since the 1950s (precipitation from the 1950s repeated during the next decades, labeled P50s) yielded a mean ELA for 2000-2009 ∼100 m lower than the actual modeled ELA, while the same ELA for unchanged temperature (labeled T50s) was only ∼40 m lower. This indicates that precipitation decrease has contributed for about 70% to the rise in ELA from the 1950s to date, while temperature has had a much smaller contribution. A slight underestimation of the modeled ELA compared to the measured mean ELA for the period 2000-2009 should be noted anew, although standard deviation of the measured ELA is smaller than the ones for modeled values, indicating that the model could be slightly too sensitive compared to measurements. A similar analysis was carried out on mass balance measurements (obtained from the difference between two digital elevation models) and estimates using the same PDD. The PDD forced with measured values from Port-aux-Français indicates an equilibrium in 1960 [START_REF] Frenot | Azorella selago Hook. Used to estimate glacier fluctuations and climatic history in the Kerguelen Islands over the last two centuries[END_REF] and a similar modeled mass balance between 1963 and 2000 (-1.13 ± 0.13 m w.e. a -1 ) as the one measured by [START_REF] Berthier | Ice wastage on the Kerguelen Islands (49 • S, 69 • E) between 1963 and 2006[END_REF] (-1.33 ± 0.90 m w.e. a -1 ). When comparing these values with the ones modeled for 2000-2009 using measured values and hypothetical cases P50s and T50s, we find that similar to ELA, 77% of the differences in surface mass balance can be explained by the decrease in precipitation.

Journal of Geophysical

We are aware that the analysis of such hypothetical series is rather simplistic and that further analysis should be needed. For instance, changes in precipitation occurrence were likely associated to variations in other climatic settings which may have had significant impact on glaciers, as, for instance, cloudiness. We analyzed the variations of this variable, which may have changed significantly during the last 50 years. Figures 8c and8d show incident shortwave and incoming longwave radiation trends that differ depending on the reanalysis. Recent reanalyses show a slight increase in incident SW radiation and a slight decrease in incoming LW radiation in recent decades. This indicates a decrease in cloud cover which certainly played a role in the

Discussion

The PDD model reproduced current and past SMB (in the early 1970s) reasonably well, showing that temperature and precipitation were the first-order variables that explained the SMB variations. PDD results indicated that the mean snow limit was about 100 m lower in the early 1970s than today. Nevertheless, an acceleration in point ablation up to high elevations and a resulting decrease in SMB gradient with elevation have been observed, which are not fully reproduced by the PDD model. Our rough estimates of SEB based on MODIS albedo showed that they result from feedback including albedo and spatial reduction in snow corniches due to lower amounts of solid precipitation. Including wind drift and snow accumulation processes is necessary for the study of fine-scale (about 10-100 m) corniche formation. However, at larger scales (1-10 km), a distributed surface heat budget would be enough to model most processes and feedback that occur at the ice cap scale.

Modeling the point SMB for 2011 using albedo estimates from MODIS imagery showed that the LW radiation budget represents 4% to 13% of the available energy for melting, while the SW radiation budget represents 58% to 64%. Consequently, turbulent heat fluxes contribute 29% to 32% to ablation. This is lower than on Patagonian glaciers [i.e., [START_REF] Schneider | Energy balance in the ablation zone during the summer season at the Gran Campo Nevado Ice Cap in the Southern Andes[END_REF] and less than estimates based on summer measurements on Ivory Glacier, New Zealand (46%) [START_REF] Hay | The synoptic climatology of ablation on a New Zealand glacier[END_REF], and on Ampere Glacier (42%) [Poggi, 1977a]. Similar values of turbulent heat fluxes contribution (31%) to melt were, however, reported for the Main Divide of the Southern Alps of New Zealand [START_REF] Neale | Energy balance and synoptic climatology of a melting snowpack in the Southern Alps, New Zealand[END_REF]. These estimates were confirmed by a 1 year (2011) SEB model on Ampere Glacier, yielding mean contributions of 68% for SW, 1% for LW, and 31% for turbulent heat fluxes, although this annual modeling should be validated in the future. Albedo thus plays an important role in ablation through the SW radiation budget, which motivated us to undertake a thorough study of changes in the snow line using MODIS data and modeling.

MODIS imagery also revealed a change in the location of the snow line between 2000 and 2012. The current snow line is above 800 m asl (Figure 4) but slightly lower on the ice cap's western flank than on its eastern flank (Figure 6), due to dominant eastward moisture transport [START_REF] Garreaud | Present-day South American climate[END_REF]. MODIS also showed a decrease in the accumulation area ratio since 2000 (Figures 5 and6), despite the larger accumulation area and lower snow line altitude values in March-April 2008 and 2009 compared to other years that were probably due to the occurrence of large precipitation amounts and relatively low temperatures in 2008 and 2009 (Figure 8). As the ice cap is very flat between ∼600 m asl and its summit (1030 m asl), even a small rise in the snow line has dramatic consequences for the reduction of the accumulation area (Figure 5). The extent of the current accumulation area in a strict sense is small, and negative SMB can probably be found at high altitudes.

In fact, comparison with field data from the 1970s indicates that the snow line has risen more than 150 m in 40 years. Moreover, the SMB gradient with elevation is lower today than it was in the early 1970s. Indeed, since snow cover is almost continuously absent up to high elevations, the albedo at stakes with distinct elevations is more uniform (as it remains in the albedo range of ice almost all through the entire ablation season), leading to a lower SMB gradient with elevation.

Most of Cook Ice Cap's wastage during the last four decades is due to reduced accumulation and subsequent rise in the snow line. In comparison, glaciers of New Zealand's Southern Alps lost about 8.4 km 3 between 1976 and 2008, more than 70% of which can be explained by thinning and calving of the largest glaciers [START_REF] Chinn | Annual ice volume changes 1976-2008 for the New Zealand Southern Alps[END_REF]. Snow line increase of small and medium glaciers, measured using oblique aerial photography [START_REF] Clare | Interannual variation in end-of-summer snowlines of the Southern Alps of New Zealand, and relationships with Southern Hemisphere atmospheric circulation and sea surface temperature patterns[END_REF] and ASTER and MODIS images [START_REF] Mathieu | Detecting the equilibrium'-line altitudes of New Zealand glaciers using ASTER satellite images[END_REF][START_REF] Sirguey | Subpixel monitoring of the seasonal snow cover with MODIS at 250 m spatial resolution in the Southern Alps of New Zealand: Methodology and accuracy assessment[END_REF][START_REF] Gjermundsen | Assessment of multispectral glacier mapping methods and derivation of glacier area changes, 1978-2002, in the central Southern Alps, New Zealand, from ASTER satellite data, field survey and existing inventory data[END_REF], is responsible for less than 30% of glacial loss over the same period [START_REF] Chinn | Annual ice volume changes 1976-2008 for the New Zealand Southern Alps[END_REF]. In Patagonia, wastage of tidewater and lake-calving glaciers is also very important [START_REF] Warren | The calving glaciers of southern South America, Global Planet[END_REF]Willis et al., 2012aWillis et al., , 2012b;;[START_REF] Melkonian | Satellite-derived volume loss rates and glacier speeds for the Cordillera Darwin Icefield, Chile[END_REF].

Conclusions

To sum up, the snow line has risen with time, the accumulation area is now reduced, and the mean albedo of the ice cap is lower. These processes led to a retreat of the ice cap. Climatic forcings that explain this Journal of Geophysical Research: Earth Surface 10.1002/2014JF003329 retreat were investigated using meteorological measurements and reanalysis data (Figure 8) and showed that changes in precipitation represent the initial first-order process that explains the marked increase in ablation and decrease in accumulation since the 1970s. An analysis of the ELA (Table 4) revealed that about 70% of the rise in ELA since the 1950s was due to the decrease in precipitation. In addition, a similar analysis, but concerning mass balance values, was carried out, also indicating a significant contribution from precipitation (77%) compared to temperature.

The decrease in precipitation was due to a reduction in the frequency of precipitation events (Table 3). This was probably linked to a southward shift of the storm tracks due to more frequent positive phases of the Southern Annual Mode in recent decades, which has had important impacts on the environment [START_REF] Villalba | Unusual Southern Hemisphere tree growth patterns induced by changes in the Southern Annular Mode[END_REF][START_REF] Manatsa | Link between Antarctic ozone depletion and summer warming over southern Africa[END_REF][START_REF] Delworth | Regional rainfall decline in Australia attributed to anthropogenic greenhouse gases and ozone levels[END_REF]. Combined with a more moderate increase in temperature, this resulted in a rise in the 0 • C level, leading to a decrease in the occurrence of solid precipitation at low to intermediate elevations. The latter caused an albedo feedback, which, combined with the slight increase in SW radiation (+2% to +4% according to data from ERA-Interim and NCEP 2 reanalysis, even though these trends must be considered with caution since long-term measurements are lacking to confirm this hypothesis), led to more ablation and a rapid rise in the snow line. The subsequent decrease in elevation caused additional positive feedback on ablation due to the local increase in temperature. However, different behaviors were observed between two zones of the glacier (zone 1 and zone 2; see Figure 1b), which are visible as large snow corniches or dunes in the lee of crests. As can be seen in Figure 1b, the spatial distribution of these zones is not related to elevation but to the distance to the crests. Because rough estimates of SMB using MODIS albedo suggest that large negative SMB should be observed, it is very likely that these zones only exist because of overaccumulation. The limited extent of these local accumulation zones is visible on the MODIS images as well as in the field. These different zones were not fully reproduced by a simple approach based on changes in fluxes with elevation, and overaccumulation processes must thereby be accounted for.

As a consequence, a full modeling approach (including full SEB modeling and wind scouring) would be useful. To this end, recently initiated regional scale simulations of mass balance processes over the Kerguelen Islands with a downscaling scheme Surface Mass balance High-resolution downscaLing (SMHiL) [START_REF] Agosta | Evolution du bilan de masse de surface Antarctique par régionalisation physique et conséquences sur les variations du niveau des mers[END_REF][START_REF] Agosta | High-resolution modelling of the Antarctic surface mass balance, application to the twentieth, twenty first and twenty second centuries[END_REF]] and the regional climate model MAR [START_REF] Gallée | Development of a three-dimensional meso-𝛾 primitive equation model: Katabatic winds simulation in the area of Terra Nova Bay, Antarctica[END_REF] should provide new insights into changes in its ice bodies. First of all, we would like to thank the Editor, Associate Editor, and reviewers for their precious advice during the review of this article. We also thank Etienne Berthier, Antoine Rabatel, Julie Gardelle, and Yves Arnaud for their help with satellite imagery and the Météo France team at Port-aux-Français for providing meteorological data and information about the meteorological station. Special thanks to the IPEV logistics team and all the people who helped us during the field campaigns without whose help this study would not have been possible. We also thank the French Austral and Antarctic Territories (TAAF). NCEP 1 and 2 reanalysis data were provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their Web site at http://www.esrl.noaa. gov/psd/. ERA-40 and ERA-Interim reanalysis data were downloaded from the ECMWF data portal at http:// data-portal.ecmwf.int/. MODIS Level 1B Swath images were retrieved from the LAADS at http://ladsweb.nascom. nasa.gov/ and the Landsat images from the USGS at http://glovis.usgs. gov/. The ASTER image is courtesy of the GLIMS project.

Figure 1 .

 1 Figure 1. (a) MODIS©image of the Kerguelen Islands showing the main study sites (in red) and the location of Port-aux-Français scientific station. The location of the Kerguelen Islands is indicated by the red dot on the world map in the top right corner. (b) Field network setup in 2010 on and around Ampere Glacier. The positions of the ablation stakes in Vallon's 1970 study are also shown. Position of the accumulation point from section 3.4 is also indicated (the latter does not correspond to a stake measurement). The background image is an ASTER image acquired on 23 April 2009. Coordinates are in universal time meridian (UTM) (Zone 42S).

Figure 2 .

 2 Figure 2. Comparison between (a) a MODIS-derived albedo map and (b) an ASTER image of the Cook Ice Cap, both produced on 23 April 2009. In Figures 2a and 2b, the dotted black line represents an albedo threshold of 0.55 and the continuous black line a threshold of 0.65. The dotted blue line shows the limits of Ampere glacier. (c) Detail of the NW part of the ice cap showing the position of the snow line for a threshold of 0.55 (dotted black line) and 0.65 (continuous black line). Clouds (Cl) are in black in Figure 2a. In Figure 2b, some clouds and their shadows across the surface are visible in the center of the ice cap. Illumination of the ASTER image was corrected to avoid saturation. Coordinates are in UTM (Zone 42S).

  4. Here we used data from 21 December 2010 to 4 January 2011 and from 14 December 2011 to 30 December 2011. During this period, an AWS was placed on the Journal of Geophysical Research: Earth Surface 10.1002/2014JF003329

Figure 3 .

 3 Figure 3. Computed and measured surface height in the ablation zone during the summer field campaigns in (top) December 2010 and (bottom) December 2011. Thick lines are the computed ice (in black) and snow (in gray) surface height computed with the SEB approach.Thin black and gray lines with squares are the ice and snow surface height (respectively) computed with the PDD approach. Black and gray circles are the mean measured ice and snow surface height (respectively) at a stake farm (four stakes) located around the Ampere AWS on the ablation zone (280 m asl; see Figure1b). The red curve is the 30 min surface elevation measured with an acoustic gauge at the same location. Error bars are ±20 mm, which represent the mean standard deviation between stake measurements.

Figure 4 .

 4 Figure 4. 2011 Ampere Glacier SMB values (m water equivalent, diamonds) at stake locations (see Figure1b) plotted against surface elevation (m above sea level). Values obtained between 1970 and 1974 byVallon [1977aVallon [ , 1977b[START_REF] Vallon | Glaciologie à Kerguelen[END_REF] are also given for the purpose of comparison. Lines represent ice cap SMB derived from our PDD model, for the early 1970s (red line, average of 1970-1974) and 2011 (blue lines). The PDD SMB for 2011 is calculated from January to January (light blue line) and from May to May (dark blue line). The dotted lines represent the version of the model with no initial firn cover, while the solid lines represent the version with an initial firn cover. Uncertainty intervals (±10 m w.e.) in the initial firn accumulation are also given.

Figure 5 .

 5 Figure 5. (a) Comparison between Ampere Glacier snow line computed with the PDD model and discrete snow line elevations for the Ampere Glacier obtained from albedo maps derived from MODIS images. Error bars are given for a snow albedo limit between 0.55 and 0.65 on the MODIS images. (b) Changes in the ratio of the accumulation area (snow-covered area) to the surface area of the entire ice cap between January 2000 and August 2012, for three periods: November-December (blue), March-April (red), and July-August (green) derived from clear-sky MODIS images. Error bars are also indicated for a snow albedo limit between 0.55 and 0.65 on the MODIS images. (c) Hypsometry of the Cook Ice Cap in 2009, i.e., the percentage of the surface area of ice cap below a given elevation.

Figure 6 .Figure 7 .

 67 Figure 6. Comparison between albedo maps for (left) 9 March 2001 and (right) 26 April 2012. For key to pixel colors, see Figure 2a. The snow line is represented by black lines (dotted line for the limit between ice and fragmented snow cover and continuous line for the limit between fragmented snow cover and full fresh snow cover), and contour lines in 2000 (SRTM DEM) and 2009 (SPOT DEM) respectively are in burgundy (numbers indicate contour line elevations in m asl). The pink arrows in the 2012 image point to sites where ice tongues or termini disappeared between 2001 and 2012. Coordinates are in UTM (Zone 42S).

Figure 8 .

 8 Figure 8. Meteorological (Météo France measurements at Port-aux-Français) and climatic (NCEP Reanalysis 1 and 2, ERA-40, and ERA-Interim reanalysis) data for the Kerguelen Islands since 1948: (a) temperature 2 m above the surface, (b) precipitation, (c) shortwave (SW) radiation flux, and (d) longwave (LW) radiation flux. Temperature and precipitation trends for different periods (delimited by vertical dashed lines) are also shown.

  Field measurements and research were funded by IPEV program 1048 (GLACIOCLIM-KESAACO) and by INSU program LEFE-KCRuMBLE.

Table 1 .

 1 Characteristics of Sensors Installed on the Mortadelle and Ampere AWS: Name, Data, Accuracy (According to the Manufacturer), and Measurement Height (La Mortadelle/Ampere)

	Name of Sensor	Data	Accuracy	Height
	Young 05103 anemometer	Wind speed	±0.3 m s -1	5 m/2.1 m
		Wind direction	±3 •
	Vaisala HMP45C	Air temperature	±0.4 • C	5 m/2.1 m
		Relative humidity	±3%
	Kipp and Zonen CG3 radiometer	Incoming LW radiation	±10% a	2 m/-
	Kipp and Zonen CMP3 radiometer	Incident and reflected SW radiation	±10% a	-/2 m
	P. MECANIQUE 3029 pluviometer	Amount of precipitation	±0.5 mm	0.85 m/-
	Campbell Sc. SC100 pressure gauge	Pressure	±1.5 mbar	2 m/-
	Acoustic gauge Campbell Sc. SR50A	Surface height	±1 cm	-/1.6 m

a Accuracy for daily totals.

Table 2 .

 2 Description of Fine-Resolution ASTER and Landsat Images

	Satellite	Date	Product Reference	Resolution
	ASTER	23 April 2009	SC:AST_L1A.003:2072780070	15 m
	Landsat	11 January 2001	LE71390942001011SGS00	30 m
	Landsat	27 November 2001	LE71390942001331SGS00	30 m
	Landsat	11 March 2005	LE71390942005070ASN00	30 m
	Landsat	10 August 2008	LE71390942008223SGS00	30 m
	Landsat	24 February 2011	LE71390942011055ASN00	30 m

Table 3 .

 3 Characteristics of Precipitation Events From Météo France Measurements atPort-aux-Français a Total annual amount of precipitation, total number of precipitation events, and number of events with less than 20 mm d -1 are shown for 10 year periods since 1953.

		Mean Annual Total Number of	Precipitation	Ratio of Events
		Precipitation	Precipitation	Events	< 20 mm d -1 Compared
		(mm a -1 )	E v e n t s	< 20 mm d -1	to Total Number of Events (%)
	1953-1962	1055	2402	2298	95.7
	1963-1972	567	2113	2083	98.6
	1973-1982	743	2102	2051	97.6
	1983-1992	789	2013	1947	96.7
	1993-2002	649	1928	1888	97.9
	2003-2012	626	1814	1764	97.2

a

Table 4

 4 compares mean ELA values from the PDD model and from MODIS images during the 1950s and the 2000s. The mean modeled ELA for the period 2000-2009 was 127 m higher than the mean ELA for the period 1958-1963 (last period of glacier equilibrium

Table 4 .

 4 Comparison of Mean ELA Values Calculated by the PDD Model and Estimated From MODIS Images Over the Periods 1958-1963 and 2000-2009 a PDD values are calculated with measured precipitation and temperature from Port-aux-Français (measured), with normal precipitation but temperature from the 1950s repeated until the 2000s (T50s) and with normal temperature but precipitation from the 1950s repeated every decade until the 2000s (P50s). Uncertainty estimates for 2000-2009 mean ELA correspond to the standard deviations of centered snow line values.

	Climate Forcing	Measured	T50s	P50s
	Time period	1958-1963	2000-2009	2000-2009	2000-2009
	PDD ELA (m)	549	677 ± 182	638 ± 196	578 ± 175
	MODIS ELA (m)		745 ± 147

a

of Geophysical Research: Earth Surface 10

  .1002/2014JF003329 ice cap wastage, although this result should be interpreted with caution because (1) previous generation reanalyses displayed opposite trends and (2) there are no available long-term field measurements in this area to enable us to confirm this result.

©2015. American Geophysical Union. All Rights Reserved.

Acknowledgments