Evolution of Ossoue Glacier (French Pyrenees) since the end of the Little Ice Age

Abstract

Little is known about the fluctuations of the Pyrenean glaciers. In this study, we reconstructed the evolution of Ossoue Glacier (42${}^{\circ }$46${}^{\prime }$ N, 0.45 km${}^{\mathrm{2}}$), which is located in the central Pyrenees, from the Little Ice Age (LIA) onwards. To do so, length, area, thickness, and mass changes in the glacier were generated from historical data sets, topographical surveys, glaciological measurements (2001–2013), a ground penetrating radar (GPR) survey (2006), and stereoscopic satellite images (2013). The glacier has receded considerably since the end of the LIA, losing 40 % of its length and 60 % of its area. Three periods of marked ice depletion were identified: 1850–1890, 1928–1950, and 1983–2013, as well as two short periods of stabilization: 1890–1894, 1905–1913, and a longer period of slight growth: 1950–1983; these agree with other Pyrenean glacier reconstructions (Maladeta, Coronas, Taillon glaciers). Pyrenean and Alpine glaciers exhibit similar multidecadal variations during the 20th century, with a stable period detected at the end of the 1970s and periods of ice depletion during the 1940s and since the 1980s. Ossoue Glacier fluctuations generally concur with climatic data (air temperature, precipitation, North Atlantic Oscillation, Atlantic Multidecadal Oscillation). Geodetic mass balance over 1983–2013 was $-$1.04 $\pm$ 0.06 w.e.a${}^{-\mathrm{1}}$ ($-$31.3 $\pm$ 1.9 m w.e.), whereas glaciological mass balance was $-$1.45 $\pm$ 0.85 m w.e. a${}^{-\mathrm{1}}$ ($-$17.3 $\pm$ 2.9 m w.e.) over 2001–2013, resulting in a doubling of the ablation rate in the last decade. In 2013 the maximum ice thickness was 59 $\pm$ 10.3 m. Assuming that the current ablation rate remains constant, Ossoue Glacier will disappear midway through the 21st century.

Introduction

The Pyrenees are a mountain range in southwestern Europe spanning $\sim 430$ $\mathrm{km}$ from the Bay of Biscay (Atlantic Ocean) to the Mediterranean Sea. According to regional climate model projections, the thickness and duration of its snowpack could decline over the 21st century . However, analysis of snow depth observations over 1985–1999 in the Spanish Pyrenees showed contrasting trends, with increasing snow depth above 2200 $\mathrm{m}$ elevation and decreasing snow depth below 2200 $\mathrm{m}$ . Tree-ring time series from living trees and in situ relict samples, collected at elevations of 2200–2450 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, have allowed the reconstruction of 1260–2005 summer temperatures in the Pyrenees. The data confirmed warming in the twentieth century . For recorded data, the longest meteorological time series in the French Pyrenees began in 1882 at an astronomical observatory located on the Pic du Midi (2862 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). A mean annual temperature increase of 0.83 ${}^{\circ }$C was observed over 1882–1970 with a significant decrease in the mean annual diurnal temperature range (2.89 ${}^{\circ }$C per century), mainly due to a significant increase in the daily minimum temperature . Recent work on data homogenization within the framework of the Pyrenean Climate Change Observatory depicts a uniform warming for the massif over the last 60 years, and highlights a significant warming signal from the 1980s onwards .

The Pyrenees hosts the southernmost glaciers in Europe, all below 43${}^{\circ }$ N latitude. Their small sizes ($<\mathrm{1}$ ${\mathrm{km}}^{\mathrm{2}}$), relatively low elevations, and southern locations make them particularly vulnerable to climate warming . Pyrenean glaciers are strongly out of balance with regional climate and are quickly retreating . While Pyrenean glaciers are in jeopardy, little is known about their evolution since the end of LIA. Their comparisons with other mountain range glaciers (e.g. Alps) are rare, and hampered by fragmented data sets . Due to the paucity of meteorological measurements, especially at high altitude, Pyrenean climate proxy records are useful to complete past climate fluctuations at secular scales. A part of this is that glaciers are considered robust climate proxies ; their reconstruction may provide further independent evidence that the climate is changing. More generally, however, retreat of Pyrenean glaciers could affect local ecosystems by diminishing the beta diversity in Pyrenean streams . Furthermore, natural patrimony and the visual perception of the high mountain landscape could also be irrevocably affected .

Ossoue Glacier (42${}^{\circ }$46${}^{\prime }$ N, 0.45 ${\mathrm{km}}^{\mathrm{2}}$) is the second largest glacier in the Pyrenees. In comparison with that of other Pyrenean glaciers, the evolution of Ossoue Glacier is well documented, with observations starting at the end of the 19th century. These include historical data sets, topographic maps, aerial images, and stake measurements.

The objective of this paper is to reconstruct the evolution of Ossoue Glacier based on the available data to provide further information:

• on the evolution of a Pyrenean glacier since the end of the LIA

• on the comparison between Pyrenean and Alpine glacial evolution

• on the potential climate drivers of Ossoue Glacier

• on the likely evolution of Ossoue Glacier in the near future.

The first section gives a brief review of studies on glaciers in the Pyrenees (Sect. ). From there, and after describing the site of Ossoue Glacier (Sect. ), the data, and the methods used (Sect. ), we propose a reconstruction that includes various glacial metrics (i.e. glacier length, area, and mass balance) for Ossoue Glacier (Sects. ). The combination of these metrics allows us, for the first time, to depict a consistent evolution of the glacier since the LIA. We compare the timeline of Ossoue Glacier fluctuations with that of Pyrenean and Alpine glacial reconstructions (Sect.  ) and with the meteorological time series (temperature, precipitation). We also mention a possible connection with the North Atlantic Oscillation (NAO) and the Atlantic Multidecadal Oscillation (AMO) (Sect. ). We use ground penetrating radar (GPR) measurements collected in 2006 to estimate the ice depth in the upper part of the glacier and to estimate the evolution of the glacier thickness in the coming decades (Sect. ). The implications and the limitations of the results are discussed in section  , and summarized in section  .

Glacier studies in the Pyrenees

The last favourable period to glacier development in the Pyrenees was the Little Ice Age (LIA), which occurred between the 14th and 19th centuries . LIA climate cooling in the Pyrenees led to the formation and advancement of glaciers in 15 massifs, in which there are up to 100 cirques . In the middle of the 19th century, after the respective advance and recession phases, the Pyrenean glacier fronts reached positions close to their maximum LIA extent. At that time, the area of Pyrenean glaciers is estimated to have been slightly over 20 ${\mathrm{km}}^{\mathrm{2}}$ . Since then, their area covered 8 ${\mathrm{km}}^{\mathrm{2}}$ in 1984 , 6 ${\mathrm{km}}^{\mathrm{2}}$ in 2004 and approximately 3 ${\mathrm{km}}^{\mathrm{2}}$ in 2013 .

Due to their remote locations and small sizes, Pyrenean glaciers have not benefited from long-term glaciological studies . That said, early topographic measurements were made by “Pyreneists”, alpinists who became enthusiasts in the exploration and observation of the Pyrenees . However, it was not until the Commission internationale des glaciers (CIG) was created in 1894 in Zürich, which led thereafter to the present-day International Association of Cryospheric Sciences (IACS) , that the situation slightly improved. Its first president, the Swiss scientist François-Alphonse Forel, promoted the organized monitoring of glaciers in the Pyrenees for comparison to the evolution of the glaciers in the Alps . Prince Roland Bonaparte established and communicated to the Commission the first regular observations of glacier frontal variations between 1874 and 1895 . Next, Gaurier monitored the glaciers over the period 1904–1927, which was interrupted by World War I . On the French side of the Pyrenees, Eaux et Forêts, the French national agency in charge of forest and water management, took over the measurements in 1932 and after World War II, during the 1945–1956 period . At the end of the 1970s, under the initiative of François Valla from the Centre Technique du Génie Rural et des Eaux et Fôrets and with the support of the Parc National des Pyrénées, the first, to our knowledge, mass balance measurements in the Pyrenees were performed at Ossoue and Taillon glaciers between 1978 and 1984 (with only qualitative data taken in 1983 and 1984) . This initiative led to the creation of the Groupe d'Etudes des Glaciers des Pyrénées (GEGP), a collaborative group comprising the Institut National de l'Information Géographique et Forestière (IGN) and researchers at Pau University. Two topographic maps, dated from 1948 and 1983 (and shown below), were produced by the GEGP . However, this group lasted only a few years, so that, for the period between 1957 and 2001, only raw terrestrial and aerial images are available for reconstructing glacier front and area variations. Since 2001, a group of volunteer glaciologists called the Association Moraine have performed regular glaciological field measurements . On the Spanish side, the institutional program Evaluación de Recursos Hídricos Procedentes de Innivación (EHRIN) has monitored Spanish glaciers since the 1990s. Since 1991, this program has collected an uninterrupted glaciological mass balance time series of the Maladeta Glacier still ongoing;.

In spite of all these efforts, observations of the Pyrenean glaciers remain scarce and irregular. Hence, there are few available reconstructions of glacier evolution since the LIA, and quantitative studies are even rarer. A brief review of Pyrenean glacier evolution is given in . On the Spanish side of the Pyrenees, the ice-covered area has decreased by 74 % since the end of the LIA . In comparison, Alpine glaciers lost 35 % of their total area from 1850 until the 1970s, and almost 50 % by 2000 . Field measurements completed by early maps, paintings, and terrestrial and aerial photographs have allowed the reconstruction of the fluctuations of the Taillon, Maladeta, and Coronas glaciers throughout the 20th century . The results of these studies are consistent with general glacial recession since the LIA. Each glacier experienced alternating periods of strong recession with periods of stability or limited readvance. In particular, there seems to be a common period of strong recession after 1850, a period of readvance or stability between 1960 and 1980, and a period of strong recession from the mid-1980s until now.

One of the main drivers of these glacial changes since the LIA is regional temperature increase . Periods of low precipitation were identified without evident trends; however, research has been lacking in order to identify potential connections to larger-scale atmospheric patterns. Generally, there are correlations with the NAO affecting the snow accumulation in the Pyrenees during the second half of the 20th century, in particular at high elevations . The AMO was more recently identified as a possible driver of multidecadal variations in river flow and precipitation in southwest France, including the Pyrenees .

Local topo-climatic effects, such as avalanching, wind-drifted snow, or shading, may significantly influence accumulation and ablation processes. In the Pyrenees, these local influences are expected to have introduced spatial disparities in ice shrinkages; in particular, to have promoted steep north and northeast glacier cirques, located below the highest summits .

Study site

Ossoue Glacier (42${}^{\circ }$46${}^{\prime }$15${}^{\prime \prime }$ N, 0${}^{\circ }$08${}^{\prime }$40${}^{\prime \prime }$ W) is located in the central part of the Pyrenees, beneath the border pass of Cerbillona. It belongs to the Vignemale Massif, which owes its name to the eponymous highest peak of the French Pyrenees (3298 $\mathrm{m}$ a.s.l.) (Figs.  and ). The glacier is an east-facing cirque. Its bedrock comprises metamorphic limestone ridges and quartzite rocks from the Devonian period.

Meta-data of Ossoue Glacier topographic survey. Ci indicates the contour interval (m) of the topographic maps. The last column refers to random errors calculated for each type of metric measurements. For volumetric measurements, we give the random error at the elevation associated with the DEM.

Metric	Period	Method	Source characteristics	Institution	Estimated error
				(surveyed by)
Length	1850	Moraine observation	Glacier deposits	Association Moraine	10 $\mathrm{m}$
variations	1885	Photointerpretation	Photograph	J. Vallot	10 $\mathrm{m}$
	1889–1895	Length measurements (field)	–	R. Bonaparte	10 $\mathrm{m}$
	1904–1928	${}^{\prime \prime }$	–	L. Gaurier	10 $\mathrm{m}$
	1935–1953	${}^{\prime \prime }$	–	Eaux et Fôrets	5 $\mathrm{m}$
	1957	Photointerpretation	Aerial image	IGN	10 $\mathrm{m}$
	1962; 1970	${}^{\prime \prime }$	Photograph	J. M. Grove; J. Jolfre	5 $\mathrm{m}$
	1982–1986; 1990	${}^{\prime \prime }$	${}^{\prime \prime }$	GEGP; B. Clos	3 $\mathrm{m}$
	1995	${}^{\prime \prime }$	Aerial image	IGN	3 $\mathrm{m}$
	2001–2013	Field observation	–	Association Moraine	1 $\mathrm{m}$
Area	1850	Moraine contour digitalization	Glacier deposits	Association Moraine	10 $\mathrm{ha}$
variations	1851	Glacier contour digitalization	Etat-Major map	French army	20 $\mathrm{ha}$
	1924; 1948	${}^{\prime \prime }$	Aerial images	IGN	4 $\mathrm{ha}$
	1950; 1953	Field measurements	–	Eaux et Fôrets	5 $\mathrm{ha}$
	1983; 1988; 1992	Glacier contour digitalization	Aerial images	IGN	4 $\mathrm{ha}$
	2002–2011	Topographic survey	GPS	Association Moraine	2 $\mathrm{ha}$
Height	1881–1895	Height measurements (field)	Artificial cave	H. Russell	0.8 $\mathrm{m}$
variations	1901–1913	${}^{\prime \prime }$	${}^{\prime \prime }$	L. Gaurier	0.5 $\mathrm{m}$
at Villa	1927; 1937	${}^{\prime \prime }$	${}^{\prime \prime }$	L. Gaurier	0.5 $\mathrm{m}$
Russell	1945–1953	${}^{\prime \prime }$	${}^{\prime \prime }$	Eaux et Fôrets	0.5 $\mathrm{m}$
	1983–1987	${}^{\prime \prime }$	${}^{\prime \prime }$	GEGP	0.5 $\mathrm{m}$
	2002–2013	${}^{\prime \prime }$	${}^{\prime \prime }$	Association Moraine	0.5 $\mathrm{m}$
Volumetric	1924	Terrestrial photogrammetry (DEM)	$\mathrm{1}:\mathrm{20\hspace{0.17em}000}$ map; Ci $=$ 20 $\mathrm{m}$	A. Meillon	8.6 $\mathrm{m}$
measurements	1948	Airborne photogrammetry (DEM)	$\mathrm{1}:2500$ map; Ci $=$ 2 $\mathrm{m}$	IGN; GEGP	2 $\mathrm{m}$
	1983	${}^{\prime \prime }$	$\mathrm{1}:2500$ map; Ci $=$ 2 $\mathrm{m}$	IGN; GEGP	1.6 $\mathrm{m}$
	2006	Topographic survey (DEM)	DGPS; base $<$ 1 $\mathrm{km}$	Sissyphe-EGID	1.5 $\mathrm{m}$
	2006	${}^{\prime \prime }$	GPR; 50 Mhz antenna	–	6 $\mathrm{m}$
	2013	${}^{\prime \prime }$	DGPS; base $<$ 40 $\mathrm{km}$	GEODE-CESBIO	0.6 $\mathrm{m}$
	2013	Satellite photogrammetry (DEM)	Pléiades stereo pair	CNES	1.8 $\mathrm{m}$

Ossoue Glacier is the largest glacier of the French Pyrenees and had an area of 0.45 ${\mathrm{km}}^{\mathrm{2}}$ in 2011. It is characterized by a large plateau on the upper part (mean elevation 3105 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, elevation range 3030–3200 $\mathrm{m}$). The plateau constitutes two-thirds of the overall area, and is located on a gentle slope (8${}^{\circ }$), while the lower part of the glacier has a steeper slope ($>20$${}^{\circ }$). Therefore, the elevation distribution along the 455 $\mathrm{m}$ elevation range is characterized by a relatively high median value (3076 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ in 2013). Ossoue Glacier has typical “alpine morphology”, being significantly longer (1400 $\mathrm{m}$) than it is wide (400 $\mathrm{m}$), and terminating in a double tongue.

Top: distribution of the present-day Pyrenean glaciers by mountain massifs: 1. Balaïtous; 2. Infierno; 3. Vignemale; 4. Gavarnie Monte Perdido; 5. Munia; 6. Posets; 7. Perdiguère; 8. Aneto; 9. Mont Valier. Bottom: Vignemale Glacieret: 1. Clot de la Hount. Glaciers: 2. Oulettes de Gaube; 3. Petit Vignemale and Ossoue. Vanished glaciers: 4. Spanish Montferrat; 5. Tapou; 6. French Montferrat. We note that the vanished glaciers were oriented to the southwest and east. Clot de la Hount is northwest-oriented and its area is less than 0.01 ${\mathrm{km}}^{\mathrm{2}}$ (2011). North-oriented glaciers Oulettes de Gaube, 0.13 ${\mathrm{km}}^{\mathrm{2}}$ (2011), and Petit Vignemale, 0.03 ${\mathrm{km}}^{\mathrm{2}}$ (2011), were one unique glacier until 1888 . Coordinate system: UTM 31${}^{\circ }$ N.

[Figure omitted. See PDF]

Ossoue Glacier is 150 $\mathrm{km}$ from the Atlantic coast and is thus under the influence of the North Atlantic westerlies, which bring abundant precipitation . The closest meteorological station (Gavarnie, 11 $\mathrm{km}$, 1380 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) recorded a mean annual temperature of 7.68 ${\mathrm{C}}^{\circ }$ and a mean precipitation of 1450 $\mathrm{mm}$ over 1992–2012. The glacier is fed mainly during winter by direct precipitation and wind-blown snow. Avalanching is most likely not a significant source of nourishment for Ossoue Glacier. The surrounding crest walls exhibit limited surfaces propitious for snow interception. Thus, Ossoue Glacier carries little debris on its surface, and topographic shading is quite limited. Dust particles are frequently observed on the snow surface, which likely affects glacier albedo and snowmelt in summer. The first day of the local hydrological year is fixed on 1 October . The melting period generally extends from the end of May to the beginning of October. We thus consider the hydrological summer during JJAS. Moulins are often observed during that period in the glacier upper area.

Ossoue Glacier volumetric variations ($\mathrm{\Delta }{V}_{\text{ice}}$ in ${\mathrm{km}}^{\mathrm{3}}$) and associated cumulative geodetic mass balance $B_{\text{geod.PoR}}$ (in $\mathrm{m}\mathrm{w}.\mathrm{e}.$) and geodetic mass balance rate $B_{\text{geod.a}}$ (in $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$). We considered a mean density of $d=850\pm \mathrm{\hspace{0.25em}50}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ for 1948–1983 and 1983–2001, and $900$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ otherwise. The term $\stackrel{\mathrm{‾}}{{\mathit{\sigma }}_{\text{geod.total.a}}}$ refers to the annualized random error.

Period of records	1924–1948	1948–1983	1983–2013
(PoR)	(24)	(35)	(30)
$\mathrm{\Delta }{V}_{\text{ice}}$	$-$0.0324 ${\mathrm{km}}^{\mathrm{3}}$	$+$0.0044 ${\mathrm{km}}^{\mathrm{3}}$	$-$0.0219 ${\mathrm{km}}^{\mathrm{3}}$
$B_{\text{geod.PoR}}$	$-$34.1 $\mathrm{m}\mathrm{w}.\mathrm{e}.$	$+$4.8 $\mathrm{m}\mathrm{w}.\mathrm{e}.$	$-$31.3 $\mathrm{m}\mathrm{w}.\mathrm{e}.$
$B_{\text{geod.a}}$	$-$ 1.42 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$	$+$ 0.13 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$	$-$ 1.04 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$
$\stackrel{\mathrm{‾}}{{\mathit{\sigma }}_{\text{geod.total.a}}}$	$\pm \mathrm{\hspace{0.25em}0.37}\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$	$\pm \mathrm{\hspace{0.25em}0.07}\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$	$\pm \mathrm{\hspace{0.25em}0.06}\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$

Ossoue Glacier mass balance time series measured by glaciological methods (in $\mathrm{m}\mathrm{w}.\mathrm{e}.$). ${\text{End}}_{\mathrm{w}}$ and ${\text{End}}_{\mathrm{s}}$ refer to the end of winter and the end of summer, respectively, in the floating-date system . $B_{\text{glac a}}$ means annual glaciological mass balance, $B_{\mathrm{w}}$ refers to the winter mass balance, and $B_{\mathrm{s}}$ refers to the summer mass balance.

Year	${\text{End}}_{\mathrm{w}}$	${\text{End}}_{\mathrm{s}}$	$B_{\mathrm{w}}$	$B_{\mathrm{s}}$	$B_{\text{glac a}}$
2002	30 May	3 Oct	2.09	$-$2.93	$-$0.85
2003	6 Jun	27 Sep	3.23	$-$4.11	$-$0.88
2004	29 May	10 Oct	3.55	$-$4.77	$-$1.22
2005	28 May	25 Sep	2.58	$-$5.07	$-$2.49
2006	25 May	8 Oct	1.95	$-$4.66	$-$2.71
2007	25 May	20 Oct	2.66	$-$4.04	$-$1.38
2008	6 Jun.	12 Oct	3.24	$-$3.35	$-$0.12
2009	30 May	12 Oct	3.15	$-$4.78	$-$1.63
2010	29 May	9 Oct	3.01	$-$3.47	$-$0.46
2011	28 May	9 Oct	2.12	$-$4.56	$-$2.44
2012	26 May	14 Oct	2.36	$-$5.78	$-$3.42
2013	7 Jun	6 Oct	3.79	$-$3.57	$+$ 0.23
Mean 2001–2013	–	–	2.81	$-$4.26	$-$1.45

Ossoue Glacier was irregularly monitored throughout the 20th century, but has been quite well monitored since 2001 (Table ).

Data sets and methods

Topographic surveys

Early sources

As is usual in glacier reconstructions, our data come from various sources (Table ). Moraines allow us to determine the glacier extent at dates estimated to be close to the end of the LIA. The testimony of Henri Passet establishes that the glacier reached the summit of the left lateral moraine in 1865 . A photograph taken in 1885 by Joseph Vallot provides evidence that Ossoue Glacier was still close to its moraines at this date. The Etat-Major map edited in 1851 by the French army also provides similar evidence. The map has an estimated accuracy of 15 $\mathrm{m}$ in planimetry. Two elevation points located on the front of the glacier are marked at 2458 and 2471 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ Currently, both points are located on the glacier moraine. At these locations, the present elevations are 2447 and 2491 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ It is remarkable that the differences in elevation are only 11 and 20 $\mathrm{m}$, which gives us further confidence in the fact that the glacier front was actually in contact with its moraine at the middle of the 19th century, i.e. at the estimated end of the LIA in the Pyrenees .

Length measurements were based on field observations reported by various authors e.g. Commission internationale des glaciers and complemented with estimations from photographs or aerial images.

The Villa Russell is a cave, accessible from the glacier at 3201 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ It was extruded by Henry Russell and his employees in 1881 (Fig. , Table ). Vertical measurements between the glacier surface and the cave threshold were made beginning in 1882.

We collected three paper topographic maps from 1924, 1948, and 1983 (Fig.  and Table ). The map dated from 1924 is a $\mathrm{1}:\mathrm{20\hspace{0.17em}000}$ scale topographic map with 20 $\mathrm{m}$ contour lines. It was created by Alphonse Meillon, a Pyrenean topograph-alpinist from the Club Alpin Français, and Etienne de Larminat, a military cartographer . Its implementation involved both field measurements and triangulation from photographs. Most of the photographs were terrestrial photographs, but, in a unique collaboration, military aerial photographs were also used to fill the information gaps . The maps from 1948 and 1983 feature 2 $\mathrm{m}$ contour lines, and were drafted by GEGP (Sect. ). Elevation contour lines were generated by manual restitution from stereoscopic airborne photographs . Both maps have a $\mathrm{1}:2500$ scale and were projected in Lambert 3 (the official French coordinate projection system until 2001).

We also collected summer aerial photographs, which date from 1924, 1948, and 1983, made available from the IGN in digital format. The latter two photographs exhibit crevasse features that match the aforementioned topographic maps, which indicates that they were the stereoscopic images used to generate the contour lines in the first place. We used these photographs to delineate the glacier outline and compute the glacier area, because we found that the glacier outline on the topographic maps was either incomplete or inaccurate. We also used the Etat-Major map (dated 1851) to compute the glacier area. We preferred the outline derived from the moraine position to that of this map in determining the glacier outline in and around the 1850s.

Orthorectification, photointerpretation, length, and area measurements (based on graphical or digital sources) were performed in a GIS (ArcGis 10.2 from Esri^©).

On the three digitized maps, contour lines were densely sampled to generate close elevation points. Two (m) DEMs were generated by interpolation, based on a discretized thin plate spline technique Anudem 5.3,.

Recent surveys

The length, area, and height (Villa Russell) measurements have continued in the 2000s.

To complete the historic DEM time series, two DGPS surveys (DGPS receivers Trimble GEO XH 2008 and 6000) were performed on 3 September 2011 and 6 October 2013 (Table ). Post-corrections, based on a 40 km distant base from the French geodetic permanent network (RGP), were applied. Two (m) DEMs were generated from the elevation point canvas, applying the same interpolation method previously mentioned. The estimated random error on the DGPS DEM is 0.6 $\mathrm{m}$.

A Pléiades stereo pair was acquired from over Ossoue Glacier on 23 September 2013. Pléiades stereo images have been successfully used to measure geodetic mass balances of mountain glaciers to within an accuracy of about 1 m . We generated a 2 $\mathrm{m}$ horizontal resolution DEM with 1 $\mathrm{m}$ vertical resolution and 1.8 $\mathrm{m}$ vertical accuracy GeoView^© software 6.6,. In this study, the Pléiades DEM was used to generate the surface elevation of the deglaciated margin between 1983 and 2013. To map the differences in surface elevation on the glacier between 1983 and 2013, the DGPS surface elevation was used because it was acquired closer to the end of the summer.

Geodetic mass balances

The DEMs generated for 1924, 1948, 1983, 2011, and 2013 allowed us to establish a geodetic mass balance over an 89 year period. Consecutive DEMs were subtracted on a pixel by pixel basis. Volume changes derived by differentiating DEMs is based on the following equation : $$\mathrm{\Delta }V=r^{\mathrm{2}}\sum _{k=\mathrm{1}}^K\mathrm{\Delta }{h}_k,$$ where $K$ is the number of pixels covering the glacier at the maximum extent, $\mathrm{\Delta }{h}_k$ is the elevation difference at pixel $k$, and $r$ is the pixel size (2 $\mathrm{m}$ in this study).

Since we have very little information on the generation of the maps based on terrestrial (1924) or aerial photogrammetry (1948, 1983), DEMs were assessed on stable terrain following the technical recommendations given in . A GCPs data set was generated from DGPS points, collected on 23 October 2013, on the frontal margin of the glacier, i.e. on a snow- and ice-free bedrock surface. DEMs were not horizontally shifted, given the good absolute localization of the sources (5 $\mathrm{m}$ for 1924, 2 $\mathrm{m}$ for 1948 and 1983), and the limited surface cover outside the glacier that would be needed to perform such an adjustment.

The annualized geodetic mass balance $B_{\text{geod.a}}$ was calculated through the following formula : $$B_{\text{geod.a}}=\frac{\mathrm{\Delta }{V}_{\text{PoR}}}{\overline{S}}.\mathit{\rho }.\frac{\mathrm{1}}{N},$$ where $\overline{S}$ is the average glacier area of the two survey dates assuming a linear change through time, and N is the number of years in the period of record (PoR). Between 1948 and 2001, we used a mean density $\mathit{\rho }$ of 850 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ with an uncertainty range of $\pm$50 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ . Before 1948, we considered a nearly absent firn zone. Since 2001, the glaciological measurements indicate that the glacier summer surface is almost exclusively bare ice. In both cases, we considered a mean density of 900 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. For details on the error estimations, please refer to the Supplement.

Photo comparison of Ossoue Glacier (Vignemale Massif): left, 1911 (L. Gaurier); right, 2011 (P. René).

[Figure omitted. See PDF]

Distribution of stakes at Ossoue Glacier. CNES^© image Pléiades MS-09-23-2013. UTM 31${}^{\circ }$ N projection.

[Figure omitted. See PDF]

Map designed by Meillon and Larminat, with focus on Vignemale glaciers, 1933 edition (glacier data are from 1924).

[Figure omitted. See PDF]

Glaciological mass balances

Since 2001, Ossoue Glacier has been monitored by systematic winter and summer mass balance measurements performed by Association Moraine. These are available on the WGMS website (Id: 2867). The direct glaciological method was used here . The protocol was similar to that used for the Saint-Sorlin and Argentière glaciers in the Alps , and followed the technical recommendations of the GLACIOCLIM observation network . At eight sites (Fig. ), the winter and annual mass balances were determined by two specific methods: (1) the end-of-winter snow depth, with respect to the previous summer surface, measured using snow probes, and near-surface snow density, calculated by drilling and weighing calibrated cores; and (2), the annual mass balance, determined by inserting 10 $\mathrm{m}$ ablation stakes (five 2 $\mathrm{m}$ sections) into the ice. Summer ablation measurements were repeated once a month until a date close to the beginning of the next hydrological year, according to the floating-date system .

These point observations were spatially integrated using an area extrapolation method. The glacier surface was divided into eight polygons centred at each ablation stake. The polygon borders were determined through empirical considerations based on field observations, elevation aspect, and mean slope (Fig. ). Further details are given in the Supplement. The winter mass balance at a specific site $k$ can be expressed as $$b_{\mathrm{w}.k}={\mathit{\rho }}_{\text{snow}.k}{h}_{\text{snow}.k},$$ where ${\mathit{\rho }}_{\text{snow}.k}$ is the density calculated at site $k$, and $h_{\text{snow}.k}$ is the snow depth accumulated during winter on the previous summer surface.

The glacier-wide winter mass balance $B_{\mathrm{w}}$ was obtained by summing the contribution from each polygon: $$B_{\mathrm{w}}=\sum _k{b}_{\mathrm{w}.k}{W}_k,$$ where $W_k$ is the fractional surface area of the polygon $k$ within the glacier. The $W_k$ values were updated in 2006 and 2011 to reflect the evolution of the glacier geometry.

The annual mass balance was calculated using the same spatial integration method. If the field operator noted the disappearance of the winter snow layer and the presence of older firn from a previous year, a density of 600 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ was applied to that layer. If ice was observed, a constant density ${\mathit{\rho }}_{\text{ice}}=900$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ was used. Lower density values were not used because of the continuous glacier shrinkage observed since the 1980s.

Meta-data of glacier variations from the literature used in this study for comparison between Ossoue and others Pyrenean and Alps glacier fluctuations.

Glacier name	ID. number on figure	Main metric	Distance from Ossoue	Publication
Taillon	1	Length variations	30 $\mathrm{km}$
Maladeta	2	Area variations	80 $\mathrm{km}$
Coronas	3		80 $\mathrm{km}$
Saint Sorlin	4	Glaciological mass balances	550 $\mathrm{km}$
Gébroulaz	5	${}^{\prime \prime }$	600 $\mathrm{km}$	${}^{\prime \prime }$
Argentiéres	6	${}^{\prime \prime }$	650 $\mathrm{km}$	${}^{\prime \prime }$
Mer de Glace	7	${}^{\prime \prime }$	650 $\mathrm{km}$	${}^{\prime \prime }$
Sarennes	8	${}^{\prime \prime }$	550 $\mathrm{km}$
Swiss glaciers (30)	9	Geodetic mass balances	$>700$ $\mathrm{km}$

These glaciological mass balance terms can be expressed in the following equation : $$B_{\text{glac.a}}=B_{\mathrm{w}}+B_{\mathrm{s}},$$ where $B_{\text{glac.a}}$, $B_{\mathrm{w}}$, and $B_{\mathrm{s}}$ designate the glacier-wide annual, winter, and summer mass balances, respectively.

The summer balance $B_{\mathrm{s}}$ was calculated as the difference between the two measured mass balance terms.

Correlation matrix (Spearman's ${\mathit{\rho }}_{\mathrm{s}}$) between the meteorological time series and Ossoue Glacier mass balances components after removal of linear trends calculated in a fixed date system,. Correlation values given in parentheses are based on the monthly mean values. Significant correlations ($p$ values $<$ 0.05) are marked with asterisks. The Gavarnie time series was not complete enough to perform the correlations between the annual mass balance and the annual mean temperature and precipitation over 2002–2011; we reported a no data value (ND) in these cases.

	$B_{\text{glac.a}.\text{1 Oct}}$	$B_{\text{s.1 Oct}}$	$B_{\mathrm{w}.\text{31 May}}$	Gavarnie	Pic du Midi	CRU	Gavarnie	Tarbes
Variables	Mass balance			Temperature			Precipitation
Period of record	2002–2013			1992–2012	1882–2013	1858–2013	1992–2012	1882–2012
$B_{\text{glac.a}.\text{1 Oct}}$	1	0.84*	0.65*	ND	$-$0.57	$-$0.66*	ND	0.74*
$B_{\text{s.1 Oct}}$		1	0.2	$-$0.76* ($-$0.8*)	$-$0.71* ($-$0.75*)	$-$0.72* ($-$0.8*)	–	–
$B_{\mathrm{w}.\text{31 May}}$			1	ND	$-$0.38	0.15	0.71*	0.72*

Location of the glaciers used for comparison with Ossoue Glacier fluctuations.

[Figure omitted. See PDF]

The stake measurements performed from 1979 to 1985 followed similar protocol as that described above, except that the glacier was divided into five longitudinal sectors . All the details on the error calculations are given in the Supplement.

Evolution of Ossoue Glacier area since the LIA. The glacier outlines are superposed on a multispectral Pléiades orthoimage taken on 23 September 2013. UTM 31${}^{\circ }$ N projection.

[Figure omitted. See PDF]

Comparison with other glaciers in the Pyrenees and the Alps

Ossoue Glacier metrics mentioned above were pieced together to detect the main phases of glacier fluctuations. Metric variations were not considered significant if they were within the estimated error range. We identified periods of glacier accumulation or stabilization, periods of ice depletion, and undetermined periods. We selected other glacier reconstructions from the literature (Table  and Fig. ). For the Pyrenees, we considered the following glacier reconstructions (Sect. ):

• Taillon (42.69${}^{\circ }$ N,$-$0. 04${}^{\circ }$ E, northeast-oriented, 2530–2800 $\mathrm{m}$, 0.08 ${\mathrm{km}}^{\mathrm{2}}$ in 2011),

• Maladeta (42.65${}^{\circ }$ N, 0.64${}^{\circ }$ E, northeast-oriented, 2870–3200 $\mathrm{m}$ altitude range, and 0.27 ${\mathrm{km}}^{\mathrm{2}}$ in 2011),

• Coronas (42.63${}^{\circ }$ N, 0.65${}^{\circ }$ E, southwest-oriented, 3100–3240 $\mathrm{m}$ and 0.02 ${\mathrm{km}}^{\mathrm{2}}$ in 2011, ice patch since the 2000s)

These reconstructions were mostly based on length and area variations, but also used other sources when they were available (e.g. a short stake measurement time series in Taillon Glacier at the end of the 1970s). In the Alps, we considered the main phases identified by:

• and updated by on four well-studied glaciers: Saint-Sorlin, Gébroulaz, Argentière, and Mer de Glace glaciers;

• on Sarennes Glacier;

• based on thirty Swiss glaciers.

Due to the morphological differences between the Alpine and Pyrenean glaciers, we considered only the Alpine glacier reconstructions which were based on mass variations (i.e. not on the glacier length or area).

Climatic data

To infer potential drivers for Ossoue Glacier variations, we used several climatic data sets, including three mean monthly air temperature data sets.

• The closest meteorological station to Ossoue Glacier was located at Gavarnie, in the same valley (1380 $\mathrm{m}$ elevation, 11 $\mathrm{km}$ east of the glacier, time series from January 1991 to May 2012).

• To cover the glacier reconstruction period, we also used the temperature time series recorded at the Pic du Midi observatory (42.93${}^{\circ }$ N, 0.14${}^{\circ }$ E; 2874 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, 30 $\mathrm{km}$ northeast of the glacier). The time series was homogenized, and gaps were filled from January 1881 to May 2011 . We used raw data from 2011 to October 2013.

• We extended the Pic du Midi temperature time series with the CRUTEM 4 data set (5 ${}^{\circ }$ $\times$ 5 ${}^{\circ }$ gridded version) over the period 1858–1890 . The mean annual difference between these data sets due to the elevation difference between the stations was removed from the CRUTEM data set ($-$13.07 ${}^{\circ }$C).

All temperatures time series correlated well (Spearman's $\mathit{\rho }>0.96$, $p$ values $<$ 0.05). Temperature time series were averaged over summer periods (JJAS), winter periods (NDJFMA), and hydrological years (beginning 1 October)

We also used two monthly precipitation data sets.

• Monthly precipitation was recorded at the Gavarnie Valley station simultaneously with the temperature (see above for station location; period of record: January 1991 to May 2012).

• Monthly precipitation was recorded over the period 1882–2013 at Tarbes-Ossun Météo-France station (43.18${}^{\circ }$ N 0.00${}^{\circ }$ E, 360 $\mathrm{m}$ elevation, 50 $\mathrm{km}$ northeast of the glacier). These data were homogenized until 2000 and used as raw data since then .

Both precipitation data sets showed a significant correlation (Spearman's $\mathit{\rho }=0.56$, $p$ value $<$ 0.05). We calculated annual and winter precipitations (NDJFMA). October and May precipitation values were not considered because the precipitation was as often liquid as solid in these months, and thus presumably did not contribute significantly to accumulation on the glacier.

We performed correlations between the meteorological data (temperature and precipitation) and Ossoue Glacier mass balance measurements over the 2001–2013 period. To reconcile the different recording periods of the meteorological and glaciological measurements, summer and annual mass balances were linearly interpolated by using a fixed date system 1 October,. The mass balances terms $B_{\mathrm{w}}$ were extrapolated to 31 May. Due to some missing data, it was not possible to calculate the mean annual temperature and precipitation for several months at Gavarnie over the period 2002–2012. Thus we did not correlate these data with the annual mass balance. Mass balances and meteorological time series were detrended using a linear least-squares fit prior to the correlation analysis.

We also considered the Atlantic Multidecadal Oscillation (AMO) and the North Atlantic Oscillation (NAO). Both indices represent fluctuations in the North Atlantic climate and have been successfully used in glacier–climate linkages .

• For the NAO we used a winter (DJFM) index based on the monthly 1850–1999 Climate Research Unit (CRU) data set, completed by Tim Osborn's 2000–2013 NAO Update . We applied a 5 year moving average filter to smooth the signal.

• For the AMO, we used the monthly index calculated from the Kaplan sea surface temperature data set over 1861–2009 .

Geophysical surveys

A ground penetrating radar (GPR) survey was performed on 30–31 August 2006 in the upper area of Ossoue Glacier. The GPR apparatus used was a PulseEkko 100 (Sensors and Software Inc.) with 50 $\mathrm{MHz}$ unshielded antennas. Three longitudinal profiles (W–E) running from the top along the slope transition of the glacier, and four transversal (N–S) profiles, were surveyed. From this data, a bedrock map was generated (see Supplement). From the 2013 glacier DEM and the bedrock map, we generated a 2013 glacier ice thickness map (Fig. ).

To estimate the ice thickness maps in the next decades, we generated ice thickness projections based on a static approach (Fig. ).

• We interpolated, on a 4 $\mathrm{m}$ resolution grid, the mean point mass balances measured at the stake locations over 2001-2013 on the current glacier outline (2011). The interpolation technique is the same used to generate the DEMs. This interpolation is consistent with the glaciological method. The mean rate after interpolation is $-$1.5 m w.e. a${}^{-\mathrm{1}}$, while the mean glaciological mass balance rate is $-$1.45 m w.e. a${}^{-\mathrm{1}}$.

• We made the assumption that this spatialized mass balance rate will remain constant in the next decades.

• Based on this mass loss rate, we calculated the ice depth on the glacier plateau at decadal intervals from 2013.

As it is based on the superficial mass balance only, this method estimating the future ice depth does not take into account effects from basal and internal mass balances, nor the dynamics of the glacier.

Results

Ossoue Glacier metric variations

Length variations

Our reconstruction of Ossoue Glacier front shows significant glacier retreat since the LIA, with intermittent stationary phases (Fig. ).

From 1850 to 1889, Ossoue Glacier front retreated by 346 $\mathrm{m}$ ($-$8.8 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$). During the following 15 years (1889–1904), the front position was quite stationary, retreating by 11 $\mathrm{m}$ between 1892 and 1893 and by only 9 $\mathrm{m}$ between 1899 and 1904. In the 1904–1905 period, however, Ossoue Glacier front retreated by 23 $\mathrm{m}$. The following periods were characterized by stability (1905–1911) and progression (1911–1927). In our data set, the glacier reached its most advanced position of the 20th century in 1927.

Area variations

The area of Ossoue Glacier at the end of the LIA, based on moraine locations, was 112.6 $\pm$ 10 $\mathrm{ha}$. The glacier area extracted from the Etat-Major map (dated near 1851) is 115 $\pm$ 20 $\mathrm{ha}$. Between the end of the LIA and 1924, the area of Ossoue Glacier decreased by 20 %. The area decreased by a further 10 % over the 1924–1948 period. During the 1948–1983 period, the front retreated by 315 $\mathrm{m}$ until 1963, and then advanced by 156 $\mathrm{m}$, although the changes in area over this period were low ($-$3 %). Over 1983–2002, the area decreased by 17 % with a notable width reduction on the slope transition. In the early 2000s, the area of Ossoue Glacier was less than 50 % of its area at the end of the LIA.

Changes in glacier geometry mainly occurred in the lower part of the glacier (Fig. ). In the upper part, the glacier shape remained almost unchanged until 1983. From 1983 to 2013, glacier width reduced dramatically at the slope transition between the plateau and the tongue of the glacier.

Elevations differences (m) on glacier (thin dashed line) and on deglaciated margins (thick dashed line) based on differences between consecutive DEMs. UTM 31${}^{\circ }$ N projection.

[Figure omitted. See PDF]

Mass variations

Since 1924, Ossoue Glacier has lost a mean of 60 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ over the current glacier area (Fig. ). The two periods of marked ice depletion, 1924–1948 and 1983–2013, were interrupted by a stable period between 1948 and 1983 (Table ).

Between 1924 and 1948, the glacier lost $-$1.42 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ ($-$34.1 $\pm$ 8.8 $\mathrm{m}$.w.e.). The ice depletion signal was strongest in the central part of the glacier (Fig. ).

Glacier surface elevation variations in $\mathrm{m}\mathrm{w}.\mathrm{e}.$ at Ossoue stake locations from 2001 to 2013. For details of stake locations on the glacier, see Fig.  and the Supplement. Maladeta Glacier is indicated by the dashed grey line.

[Figure omitted. See PDF]

The 1948–1983 period is the only period with observed positive geodetic mass balance variation, with a rate of $+$0.13 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ ($+$4.8 $\pm$ 2.6 $\mathrm{m}\mathrm{w}.\mathrm{e}.$). However, a notable depletion was observed on the tongue (Fig. ). The glacier advanced over a very small area (1 $\mathrm{ha}$) with a mean ice growth of 6.5 $\mathrm{m}\mathrm{w}.\mathrm{e}.$. An area of higher accumulation is localized on the lower part of the glacier, below the slope transition (Fig. ). At the end of that period, using ablation stakes, François Valla and Henri Pont measured mass gains of $+$0.81 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ in 1978; $+$0.26 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ in 1979; +0.17 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ in 1980; and 0 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ in 1981 and 1982. In 1983 and 1984, they considered the mass balances to be “negative” but did not provide quantitative information .

Over 1983–2013, the glacier lost $-$1.04 $\pm$ 0.06 w.e.a${}^{-\mathrm{1}}$ ($-$31.3 $\pm$ 1.9 m w.e.) . A marked pattern of ice depletion occurs along a longitudinal profile on the upper part of the glacier (Fig. ). This phenomenon increased glacier convexity in that zone, which was once named “Plateau des Neiges” in older maps.

Between 2001 and 2013, superficial mass loss on Ossoue Glacier given by the glaciological method is 17.36 $\pm$ 2.9 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ ($-$1.45 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$) (Table  and Fig. ). The strongest mass losses were registered at the lowest elevated stakes (stakes 7 and 8). The mass balance was negative every year since the stake measurements began except in 2012–2013 with a value of $+$0.23 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ . In 2008 the mass balance was only slightly negative.

Mains phases of variations

Taken together, all the metrics suggest a clear retreat of Ossoue Glacier since the end of the LIA (Fig. ). Based on the dates of the survey reported in the available documents, we identified several phases in the fluctuations of Ossoue Glacier:

• 1850–1890 (40 years, noted phase I in Fig. ). Between these dates, all the length variations are negative for all three survey dates (1874,1885,1889). Between 1882 and 1889, however, Russell noted some thickening and thinning at Villa Russell, with high interannual variability.

• 1890-1928 (38 years, phase II). The length variations were null from 1891 to 1899, and within the range of errors between 1892 and 1893 (i.e. no significant variations). Afterwards, length variations were negative until 1905. Between 1890 and 1894, Russell noted a period of stabilization at Villa Russell (II.a). Then, several irregular thinnings were observed for the following survey dates (1895, 1898, 1902, 1904), interrupted by a thickening between 1898 and 1901. From 1905 to 1911, no length variations were observed for the six survey dates. Between 1911 and 1921, the length increased, but the variation was not significant. From 1905, the glacier thickened at Villa Russell until 1913 (II.b). The glacier thinned significantly at Villa Russell between 1913 and 1927. The length variation was positive, but not significant for between 1921 and 1927.

• 1928–1950 (22 years, phase III). In spite of a large estimated error, the mass balance variation (1924-1948) was significant and markedly negative. Area variations were also negative between 1924 and 1948, and between 1948 and 1950. The glacier thinned at Villa Russell until 1950 (survey dates for significant variations: 1937, 1945, 1950). The length variations were negative for all survey dates in that stage (1928, 1935, 1950).

• 1950–1983 (33 years, phase IV). The area variations were small and not significant between 1950 and 1953, and between 1953 and 1983. Negative length variations were observed until 1962, and positive afterwards. The geodetic mass balance was slightly positive over the 1948–1983 period. At Villa Russell, the glacier was thickening at the surveyed dates (1952, 1953, 1967). In 1953, the glacier reached the threshold of the Villa Russell, and was above the threshold in 1967. Between 1967 and 1983, a slight thinning was observed. The glaciological mass balance was positive from 1978 to 1980, zero in 1981 and 1982, and negative in 1983 and 1984 (qualitative assessment for these two dates).

• 1983–2013 (30 years, phase V). From 1983, all the significant length and area variations observed were negative. At Villa Russell, the only positive variations were observed between 1987 and 1991, and 2007 and 2008. The geodetic mass balance was negative between the two dates of survey that define that phase. The annual glaciological mass balances have all been negative since 2001, except for the hydrological year 2012-2013.

From the above considerations, we defined three periods of ice depletion, when the metrics consistently indicate a negative trend: 1850–1890, 1928–1950, and 1983–2013 (noted I, III, and V in Fig. ). By the same reasoning, three periods are characterized by stability or slight growth: 1890–1904, 1905–1913, and 1950–1983 (noted II.a, II.b, and IV, respectively in Fig. ).

Length (m), area (ha) and thickness (m) at Villa Russell and mass changes (in $\mathrm{m}\mathrm{w}.\mathrm{e}.$) of Ossoue Glacier. Glaciological mass balances of Ossoue (orange) and Maladeta (blue) glaciers. The background colour indicates the interpreted trend of the period, according to the metric variations (see Sect. ). The ice depletion periods (I, III, V) are shown in light grey. The periods of accumulation or stability (II.a, II.b, IV) are shown in dark grey.

[Figure omitted. See PDF]

If we consider the 2013 glacier area as a common integration area for all the periods, the absolute value of the geodetic mass balance increased over 1924–1948 ($-$35 $\mathrm{m}\mathrm{w}.\mathrm{e}.$) and 1948–1983 ($+$6.1 $\mathrm{m}\mathrm{w}.\mathrm{e}.$).

Comparison of Ossoue Glacier fluctuations with other Pyrenean and Alps glacier reconstructions. The ice depletion periods are shown in red. The accumulation or stable periods are shown in blue. Periods considered as “undetermined” are shown in white. The meta-data relative to the glaciers fluctuations are given in Table . For the localization of the glaciers, please refer to Fig. . Roman numerals above Ossoue Glacier fluctuations refer to the main periods of variation identified (Sect. ).

[Figure omitted. See PDF]

Comparison with Pyrenean and Alpine glaciers variations

Pyrenean and Alpine glaciers exhibit similar multidecadal variations during the 20th century. The ice depletion was particularly intense in the 1940s and since the 1980s. The stable period detected at the end of the 1970s is also evident in all the glaciers at both mountain ranges (Fig. ). The Ossoue, Maladeta, and Coronas glaciers present consistent variations: ice depletion in the 1930s and 1940s, transition in the 1950s, stable period from the 1960s to the 1980s, and ice depletion since the 1980s. The Taillon Glacier variations appear more specific, but coincide with the ice depletion phase which had started in the 1980s. The Taillon Glacier and Ossoue Glacier fluctuations present further consistencies: around 1910 (II.b), in the 1930s (III), and around 1980 (IV). In the French Alps, two steady-state periods, 1907–1941 and 1954–1981, and two periods of recession, 1942–1953 and 1982–2013, were deduced from four glacier mass balance time series . These periods are in good agreement with Ossoue and Maladeta glacier variations. However, it seems that the glacier retreat phase of 1928–1950 (III) identified from Ossoue Glacier data started about a decade later in the Alps. In the Swiss Alps, Huss et al. (2010) detected two short periods of mass gain (1910s and late 1970s) and two periods of rapid mass loss (1940s and late 1980s to present). These variations are also consistent with Ossoue Glacier variations (periods III, IV, and V).

Comparison with meteorological time series

Correlations between Ossoue Glacier mass balance time series (2001–2013) indicate that the annual mass balance is mainly dependent on the summer mass balance, and that the winter mass balance has less influence (Spearman's $\mathit{\rho }=0.84$ for summer mass balance and $\mathit{\rho }=0.65$ for winter mass balance) (Table ).

The link between ablation and air temperature was verified at Ossoue Glacier, as shown by the following significant ($p$ value $<$ 0.05) correlations (Table ):

• between monthly summer ablation and monthly air temperature time series over 2002–2013 ($\mathit{\rho }=-0.8$ for Gavarnie valley station, $\mathit{\rho }=-0.75$ for Pic du Midi station, which is located farther from the glacier, and $\mathit{\rho }=-0.8$ for the regional CRU time series);

• between mean summer air temperature and summer-wide mass balance $B_{\mathrm{s}}$ (June–September) ($\mathit{\rho }=-0.76$ for Gavarnie, $\mathit{\rho }=-0.71$ for Pic du Midi and $\mathit{\rho }=-0.72$ for the CRU time series);

• between annual mass balance and mean annual temperature ($\mathit{\rho }=-0.66$ for the CRU time series).

The link between annual mass balance and mean annual temperature is weaker in the Pic du Midi time series than in the CRU time series ($\mathit{\rho }=-0.57$ for Pic du Midi). This may be due to the use of raw data in the Pic du Midi time series, starting from 2011, or to the limited period for glaciological mass balance records (annual mass balance measurements only began in 2001). However, due to the good correlation between the CRU and the Pic du Midi temperature data sets, we also considered that the mean annual temperature at Pic du Midi is linked to $B_{\text{glac.a}}$ over the longer 1890–2013 period. The elevation of the Pic du Midi station (2874 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) is close to that of Ossoue Glacier front (2755 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$); thus, we principally used this data set to identify temperature trends over historical periods.

Mean temperature (Pic du Midi, 2874 $\mathrm{m}$), mean precipitation (Tarbes, 374 $\mathrm{m}$), and correlations between mean temperature and precipitation time series and time (Spearman's $\mathit{\rho }$). Time ranges are based on the interpretation of the glacier metrics (Fig. ). Significant correlations ($p$ values $<$ 0.05) are marked with asterisks. For 1858–1890, the mean temperature is based on the CRU data set.

Time range	Period	Mean	Mean	${\mathit{\rho }}_{\mathrm{s}}$	${\mathit{\rho }}_{\mathrm{s}}$
(glacier periods)(number of years)	of the year	temperature	precipitation	temperature	precipitation
1858–2013 (T ${}^{\circ }$C) (131 y.)	hydrological (h) year	$-$1.1 ${}^{\circ }$C	1068 $\mathrm{mm}$	0.38*	$-$0.16
1882-2013 (P $\mathrm{mm}$)	winter (w)	–	556 $\mathrm{mm}$	–	0.1
	summer (s)	$+$5.3 ${}^{\circ }$C	–	0.28*	–
1858–1890 (T ${}^{\circ }$C)	h. year	$-$1.3 ${}^{\circ }$C	–	$-$0.6*	–
(I) (32 y.)	s.	$+$5.8 ${}^{\circ }$C	–	$-$0.36	–
1890–1894	h. year	$-$ 1.4 ${}^{\circ }$C	1071 $\mathrm{mm}$	–	–
(II.a) (4 y.)	w.	–	571 $\mathrm{mm}$	–	–
	s.	$+$4.6 ${}^{\circ }$C	–	–	–
1905–1913	h. year	$-$1.8 ${}^{\circ }$C	1165 $\mathrm{mm}$	–	–
(II.b) (8 y.)	w.	–	550 $\mathrm{mm}$	–	–
	s.	$+$4.1 ${}^{\circ }$C	–	–	–
1928–1949	h. year	$-$1.2 ${}^{\circ }$C	993 $\mathrm{mm}$	0.23	-0.19
(III) (21 y.)	w.	–	493 $\mathrm{mm}$	–	-0.2
	s.	$+$5.1 ${}^{\circ }$C	–	0.26	–
1950–1982	h. year	$-$1.4 ${}^{\circ }$C	1068 $\mathrm{mm}$	0.13	0.4*
(IV) (32 y.)	w.	–	586 $\mathrm{mm}$	–	0.28
	s.	$+$4.8 ${}^{\circ }$C	–	0.1	–
1983–2013	h. year	$-$0.4 ${}^{\circ }$C	1042 $\mathrm{mm}$	0.39*	$-$0.12
(V) (30 y.)	w.	–	567 $\mathrm{mm}$	–	$-$0.14
	s.	$+$6.1 ${}^{\circ }$C	–	0.43*	–

Climatic time series: mean annual temperature at Pic du Midi (beginning 1 October), annual precipitation at Tarbes, AMO mode and winter NAO (DJFM) anomalies. The background colour indicates the interpreted trend of the period, according to the metric variations (see Sect.  and Fig. ). The ice depletion periods (I, III, V) are shown in light grey. The periods of accumulation or stability (II.a, II.b, IV) are shown in dark grey.

[Figure omitted. See PDF]

Precipitation records at Gavarnie and Tarbes are significantly correlated with the winter mass balances ($\mathit{\rho }=0.71$ for Gavarnie, $\mathit{\rho }=0.72$ for Tarbes, which is located farther from the glacier, Table ). The link between annual mass balance and annual precipitation is significant in the Tarbes data set ($\mathit{\rho }=0.74$). Thus, we conclude that the Tarbes time series can be used to identify trends in precipitation that are linked with Ossoue Glacier fluctuations.

The mean annual temperature over the hydrological year (starting 1 October) for 1858–2013 is $-$1.1 ${}^{\circ }$C, and the mean summer temperature (JJAS) is 5.3 ${}^{\circ }$C. Both time series present a linear trend over the period ($\mathit{\rho }=0.38$ for annual, $\mathit{\rho }=0.28$ for summer). This correlation is stronger if we limit the period to 1882–2013 ($\mathit{\rho }=0.54$ for annual, $\mathit{\rho }=0.51$ for summer). The annual precipitation and the winter precipitation (NDJFMA) over 1882–2013 are 1068 and 556 $\mathrm{mm}$, respectively. No trend was observed in precipitation time series.

Map: bedrock depth as interpreted from GPR radargrams superposed on a 2013 XS Pléiades image. Numbers 1 to 3: interpretations of longitudinal radargram acquisitions. Numbers 4 to 7: interpretations of transverse radargram acquisitions.

[Figure omitted. See PDF]

Map: estimated ice depth at Ossoue Plateau in the next decades. These projections are based on the sum of the ice depth estimation (measured by GPR, see Fig. 12) and the cumulative mass balance over the period of projection. The latter is based on the interpolation values of the averaged mass balance at stake locations over 2001–2013 (data per year and stakes are given in Fig. 8). Spatial resolution is 4 m. East of the plateau limit (dashed line), the ice depth was unknown in 2013.

[Figure omitted. See PDF]

Analysis of temperature and precipitation trends over the four long periods stemming from the combination of glacier metrics (I, III, IV, and V, Fig. ) reveals four significant trends (Table ):

• The mean annual temperature over 1858–1890 may have continuously decreased. Over the same period, the mean summer temperature (JJAS) is 0.5 ${}^{\circ }$C higher than the mean summer temperature over 1858–2013 (5.3 ${}^{\circ }$C).

• Both short periods of glacier accumulation (II.a and II.b) present annual ($-$1.3 and $-$1.4 ${}^{\circ }$C) and summer (4.6 ${}^{\circ }$C and 4.1 ${}^{\circ }$C) temperatures lower than that of the means over 1858–2013 ($-$1.1 ${}^{\circ }$C and 5.3 ${}^{\circ }$C). The precipitation at Tarbes station over 1890–1994 and 1905–1913 is also higher than the mean over 1882–2013.

• The annual precipitation trend over 1950–1982 is positive ($\mathit{\rho }=0.4$) and its mean is equal to the mean precipitation over 1882–2013 (1068 $\mathrm{mm}$). Winter precipitation is higher than the mean recorded over 1882–2013 (586 $\mathrm{mm}$). The annual mean temperature ($-$1.4 ${}^{\circ }$C) and mean summer temperature (4.8 ${}^{\circ }$C) are lower than the means over 1858–2013.

• The last period considered (1983–2013) shows positive trends in both mean annual and mean summer temperature, with the highest registered mean temperatures ($-$0.4 ${}^{\circ }$C for annual and 6.1 ${}^{\circ }$C for summer).

Figure  provides insight into the possible linkage between the evolution of Ossoue Glacier and the regional-scale climate. The 1960–1980 period is characterized by a succession of negative phases in the NAO. This coincides with a period of relative glacier growth or stability (positive variations in various glacier lengths, areas, and mass balances). AMO warm phases occurred during 1860–1880 and 1940–1960, and cool phases during 1905–1925 and 1970–1990 . The AMO index presents some potential correlations with Ossoue Glacier variations: periods I and III of ice depletion, in regards to the AMO warm phases; and periods II.b and IV, in regards to the AMO cold phases (Fig. ).

Ice thickness maps

In 2006, the estimated mean ice thickness was 29.3 $\pm$ 6.3 $\mathrm{m}$ (max. 74.8 $\pm$ 10.2 $\mathrm{m}$), giving an estimate of 25 $\pm$ 6.5 $\mathrm{m}$ (max. 59 $\pm$ 10.3 $\mathrm{m}$) in 2013 (Fig. ). In 2011, another GRP survey indicated a maximum depth of 45 $\mathrm{m}$ and an average depth of 30 $\mathrm{m}$. Despite the discrepancies in ice thickness, both studies suggest similar bedrock morphologies. Moulins were explored over the 2004–2009 summers at the Ossoue Glacier plateau. The depth of the explored moulins ranged from 30 m to 41.5 m. Given that the bedrock was never reached according to the speleo-glaciologist, the ice thickness obtained by GPR is consistent with these depth measurements.

Over the subsequent decades, the mean and the maximum ice depths at the Ossoue Glacier plateau would rapidly decrease: 22 (mean) and 48 $\mathrm{m}$ (max.) in 2023, 17 and 38 $\mathrm{m}$ in 2033, 11 and 27 $\mathrm{m}$ in 2043, and only 3 and 7 $\mathrm{m}$ in 2053 (Fig. ).

Discussion

Using multiple data sets, we generated five independent time series of glacier metrics (length, area, thickness at Villa Russell, and geodetic and glaciological mass balances variations) to reconstruct the evolution of Ossoue Glacier since the end of the LIA (Fig. ). The metrics give a generally consistent chronology of glacier fluctuations since the LIA, although there are some discrepancies. We should bear in mind that the metrics do not directly reflect the same glaciological processes. The time series of frontal variations offers the best temporal resolution of the onset of glacier changes, but these changes are strongly influenced by ice dynamics. Glacier motion is dependent on mass variations of the upper part of the glacier, but the response time is largely unknown. Areal variations depend on ice thickness at the edges only. Thickness variations registered at Villa Russell are the result of accumulation and ablation variations at the northern periphery of the glacier only, which could be prone to snow drifting. Volumetric mass changes generated by the geodetic method are mostly the result of the surface energy budget but also include internal and basal mass variations, which remain difficult to estimate. Glaciological mass balances reflect the link between energy and mass budget properly, but can only be measured at a limited sample of points at the glacier surface. However, these metrics are all sensitive to glacier mass changes, with different time scales and response intensities. For instance, between 1948 and 1983, the mass balance was positive, yet frontal variations were negative until 1963. This can likely be explained by the delay in the response time of frontal response to glacier mass changes. In the case of Ossoue Glacier, it is note-worthy that the metrics over the study period reveal a consistent signal (Fig. ).

The evolution of Ossoue Glacier is consistent with the reconstructed evolutions of other Pyrenean glaciers. Some discrepancies might be due to the nature of the metrics used in the reconstruction, or due to the local topo-climatic influences. Considering an accuracy of $\pm$5 $\mathrm{year}$s, the study was able to identify two common stable periods (1905–1930 and 1955–1985), as well as two periods of marked ice depletion (1850–1900 and from the mid-1980s until now). The evidence of strong marked ice depletion found in this study for Ossoue Glacier between 1924 and 1948 ($-$1.42 $\mathrm{m}$ m.w.e ${\mathrm{a}}^{-\mathrm{1}}$) should be considered with caution given the high uncertainties in the altimetry restitution process. However, reconstructions of other Pyrenean glaciers over comparable periods tend to corroborate this result. Between approximately 1928 and 1957, the length of the Coronas Glacier decreased from 600 to 350 $\mathrm{m}$, while its area decreased from 19 to 8.6 $\mathrm{ha}$ and its equilibrium line altitude (ELA) increased from 3065 to 3122 $\mathrm{m}$ . During the 1935–1957 period, the Maladeta Glacier lost 15 $\mathrm{ha}$ ($-$0.68 $\mathrm{ha}{\mathrm{a}}^{-\mathrm{1}}$) and its length decreased by 80 $\mathrm{m}$. This retreat is assumed to be due to a warm anomaly detected in the second half of the 1940s . The Maladeta Glacier mass balance time series (1991–2013) is in good agreement with Ossoue glaciological mass balance time series over the 2001–2013 period (Fig. ). During the 1990s, the Maladeta mass balance values were slightly negative. If we compare the Ossoue geodetic mass balance during 1983–2013 ($-$1.04 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$) and the Ossoue glaciological mass balance during 2001–2013 ($-$1.45 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$), we can deduce that the 1983–2001 Ossoue ablation rate was approximately $-$0.76 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$. These results are consistent with the variation of the mean annual mass balance of the Maladeta Glacier: $-$0.2 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ over 1991–2001 and $-$1.03 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$, over 2001–2013 (Fig. ).

The comparison between Pyrenean and Alpine glacier fluctuations suggests that there is a common climatic driver governing glacier fluctuations in both mountain ranges. Ossoue Glacier seems to be anti-correlated with the NAO. Similar results were reported by and for glaciers in the southern Alps. In addition, variations in the AMO index also appear relatively similar to variations in the combined Ossoue Glacier metrics throughout the 20th century (Fig. ). This result is consistent with previous studies on the influence of the multidecadal internal variability of the North Atlantic circulation on the Northern Hemisphere climate e.g. and on Alps glacier fluctuations.

Variations of Ossoue Glacier metrics are in good agreement with meteorological data: periods of ice depletion are generally characterized by lower values of mean precipitation and temperatures (Table ). The evolution of Ossoue Glacier may be partially explained by observed trends, with a significant positive trend in 1950–1982 precipitation (a stable period for the glacier) and a significant constant rise in mean annual and summer temperature since 1983 (a period of depletion for the glacier). The 1850–1890 and 1983–2013 periods are marked by ice depletion, although the mean air temperature time series have opposite and significant trends. Frontal variations and mean air temperature variations, over the 1850–1890 interval, point to a shorter period of marked ice depletion, 1850–1874, with lesser depletion over the 1874–1890 period. By the same reasoning the selected 1928–1950 period of ice depletion may have been more pronounced in the 1937–1950 “subperiod” than over the 1928–1937 “subperiod”. In the Alps, the 1942–1950 period is characterized by extraordinarily high rates of mass loss .

The future evolution of Ossoue Glacier depends on climatic changes, but is also constrained by the remaining ice volume. Assuming that Ossoue Glacier mass balance follows the same trend as that recorded during 2001–2013, the glacier should disappear in 40 $\mathrm{years}$ (Fig. ). We anticipate that the glacier may split into two parts at the slope transition (Fig. ) in the near future. At this location the glacier may be particularly thin, and there is an abrupt change in the glacier slope. The lowest part may soon no longer be fed by the ice flow from the upper area and could thus rapidly disappear. This separation would drastically change the morphology of Ossoue Glacier from an active glacier to a glacieret or ice patch. Such glacier fragmentation has been regularly observed on Pyrenean glaciers, e.g. the neighbouring Petit Vignemale and Oulettes glaciers.

However, future evolution of Ossoue Glacier based on interpretation of Fig.  was made under several strong assumptions: (i) basal and internal mass balances were neglected; (ii) ice motion was neglected; and (iii) the mass loss in the future decades will occur at the same pace as during the last decade. In short, future work is necessary to better understand the effect of local topography on the spatial variability of glacier mass balance. This influence is expected to increase in the future as the glacier retreats .

The reconstruction and the future evolution of Ossoue Glacier does present large uncertainties, and the influence of the climate fluctuations on the glacier metric variations are complex; however, considering the current ablation rate, it seems doubtless that Ossoue Glacier will disappear halfway through the 21st century. Its large, markedly convex plateau (two-thirds of the present-day area) has allowed the accumulation of a significant amount of ice at high altitude (3105 $\mathrm{m}$) during favourable periods. On the contrary, its eastern orientation and low shading may have a large influence on the high rate of summer ablation (e.g. in comparison to the Maladeta Glacier). Henceforth, due to the limited interval range of the plateau (3030–3200 $\mathrm{m}$, slope 8${}^{\circ }$), any future rise of the lower limit of the glacier (2755 $\mathrm{m}$ in 2013) would drastically modify the responses of the metrics of Ossoue Glacier to future climate fluctuations.

Conclusions

Ossoue Glacier is one of the southernmost glaciers in Europe. Using a large archive of historical data sets and recent accurate observations, we generate consistent time series of various glacier metrics, such as length, area, elevation variations, and mass changes since the LIA at high temporal resolution. The dominant trend is a retreat over the 20th century, which was interrupted by two stable short periods, 1890–1894 and 1905–1913, and a longer stable period, 1950–1983.

The evolution of Ossoue Glacier is in good agreement with those of other Pyrenean glacier reconstructions (Maladeta, Coronas, Taillon glaciers), suggesting the possibility of long-term high-elevation climate reconstruction in the Pyrenees. The comparison between Pyrenean and Alpine glacial fluctuations highlights similar multidecadal variations during the 20th century. The ice depletion was particularly intense in the 1940s and since the 1980s, while a stable period detected at the end of the 1970s is also evident in all the glaciers from both mountain ranges. This result may suggest that there is a common climatic driver governing glacier fluctuations of both mountain ranges.

The time resolution of the generated metrics for Ossoue Glacier allows us to extract consistent glacial changes over various periods. These periods appear to be roughly in phase with hemispheric climate proxies, such as the North Atlantic Oscillation and the Atlantic Multidecadal Oscillation. The 1960–1980 stable period may be partially explained by anti-correlation to the NAO index. We found that the ablation rate may have doubled in the last decade, likely as a result of the recent climate warming. Ossoue Glacier fluctuations generally concur with climatic data, suggesting that Ossoue Glacier is a good regional climate proxy.

The eastern orientation and low shading of Ossoue Glacier make it particularly vulnerable to climate fluctuations, although its relatively high elevation has allowed the accumulation of a significant amount of ice at high altitudes. In 2013, the maximum ice thickness was 59 $\pm$ 10.3 m. Assuming that the current ablation rate stays constant, Ossoue Glacier will disappear midway through the 21st century.

The Supplement related to this article is available online at doi:10.5194/tc-9-1773-2015-supplement.

Acknowledgements

This manuscript was greatly improved thanks to the constructive comments of Chris R. Stokes (Editor), J. Ignacio López-Moreno (referee), and an anonymous referee. The authors warmly thank all the volunteers who provided great help during fieldwork, especially the members of the Association Moraine (including Florian Pinchon, Antoine Simmonet, and Sabine Ayrinhac). We acknowledge all our colleagues who helped collect the GNSS data in the field (including Stephane Binet, Frédéric Blanc, Bruno Calvino, and Vincent Cabot) and acquire the Pléiades images (including Claire Tinel and Steven Hosford) within the framework of the CNES Pléiades thematic commissioning phase and the Airbus Defence and Space “Pléiades User Group”. This work was supported by the Fondation Eau, Neige et Glace through the project CLIM Ex-PYR (www.fondation-eng.org) and by the Région Midi-Pyrénées and the University of Toulouse through the CRYOPYR project. We also thank Patrick Wagnon, who kindly sponsored the CLIM Ex-PYR project, Etienne Berthier, and Christophe Kinnard for their useful comments that improved the manuscript. Edited by: C. R. Stokes