NUMERICAL SIMULATIONS OF RECENT AND FUTURE EVOLUTION OF MONTE PERDIDO GLACIER

. Glaciers are globally retreating due to climate change, and the Pyrenees Mountain range is no exception. This study uses the Open Global Glacier Model (OGGM) to explore the dynamics of the Monte Perdido glacier, one of the largest remaining glaciers in the Pyrenees. We explored three calibration approaches to assess their performances when reproducing observed volume decreases. The first approach involved mass balance calibration using terrestrial laser scanning data from 2011 to 2022 and climate data from a nearby weather station. The second approach used terrestrial laser scanning calibration with default climate data provided by OGGM (GSWP3-W5E5). The third approach used default geodetic mass balance calibration and default climate data. By comparing these calibration strategies and analysing historical data (terrestrial laser scanning and ground penetrating radar), we obtain insights of the applicability of OGGM to this small, mild conditions, Pyrenean glacier. The first calibration approach is identified as the most effective, emphasising the importance of selecting appropriate climate data and calibration methods. Additionally, we conducted future volume projections using an ensemble of General Circulation Models (GCMs) under the RCP2.6 and RCP8.5 scenarios. The results indicate a potential decrease in total ice volume ranging from 91.60% to 95.16% by 2100, depending on the scenario. Overall, this study contributes to the understanding of the Monte Perdido glacier’s behaviour and its response to climate change through the calibration of the OGGM, while also providing the first estimate of its future melting under different emission scenarios.


Introduction
Glaciers are highly sensitive indicators of recent climate variations (Beniston, 2003;Grunewald and Scheithauer, 2010).Current assessments of the Intergovernmental Panel on Climate Change (IPCC) have highlighted that changes in temperature and precipitation have resulted in global glacier retreat since the 1950s that is unprecedented in the last 2000 years (IPCC, 2021).
Glaciers in the Pyrenees are currently in a critical situation, with clear evidence of very advanced stages of degradation (Rico et al., 2017;Vidaller et al., 2021).Due to their small dimensions, glaciers in the Pyrenees have minimal impact on water resources and global albedo feedback (López- Moreno et al., 2020).However, they hold scientific and touristic value while carrying strong cultural heritage (García-López et al., 2021;Moreno et al., 2021;Serrano Cañadas, 2023).Therefore, their melting represents a significant event, symbolising the wider consequences of climate change.
Glaciers, characterized by compact, perennial ice, experience mass gain through snow accumulation and mass loss during ablation, primarily through surface melting (van der Veen, 2013;Eis, 2020).The balance between accumulation and ablation determines a glacier mass fluctuations, with retreat occurring when ablation surpasses accumulation.
Glacier dynamics of mass balance respond to climate fluctuations on longer time scales rather than immediately (Huston et al., 2021).Consequently, the advance or retreat of a glacier is not only determined by the weather of a single year but is a response to cumulative forcings from many years (Huybers and Roe, 2009).Furthermore, there are other processes that influence glacier evolution such as avalanches, being sheltered from dominant winds, debris cover thickness, slope of the ice surface, or rocky outcrops that may appear and enhance incoming long-wave radiation (López- Moreno et al., 2019).

PAPER ACCEPTED. PREPRINT VERSION
Given the urgency of climate change and glacier retreat, there is significant motivation to study glacier dynamics through a modelling approach, allowing for predictions of future volume trends given climatic and geographic inputs.Hence, this study aims to explore the performance of a glacier model and its practical implementation for the Monte Perdido glacier, one of the largest remaining glaciers in the Pyrenees (Vidaller et al., 2021), with noticeable thinning observed in recent years (López-Moreno et al., 2019).

Study area
The Monte Perdido glacier, located in the Ordesa and Monte Perdido National Park in the Central Spanish Pyrenees (42.6806°N, 0.0375°E), consisted of two ice bodies until 2021: the upper and lower glaciers (Fig. 1).Both bodies are north facing and lie beneath the Monte Perdido Peak (3355 m a.s.l.).The mean elevations of the upper and lower ice bodies are between 3110 and 2885 m a.s.l., respectively (Julián and Chueca, 2007).In 2022, the lower Monte Perdido glacier experienced a division, resulting in the formation of two separate ice bodies (Fig. 1).PAPER ACCEPTED.PREPRINT VERSION

Data analysis tools
To simulate and analyse the Monte Perdido glacier, we utilised the Open Global Glacier Model (OGGM), which is a Python based open-source model (Maussion et al., 2019).We employed version 1.6.0 of OGGM, which was released on March 10, 2023 (Maussion et al., 2023).With the glacier outlines, topographical data, and climate data at a reasonable resolution, the model can estimate the total ice volume of the glacier and simulate its dynamic evolution in response to different climate forcings (Maussion et al., 2019).To compute the ice thickness, the model uses an ice thickness inversion method based on Farinotti et al., (2009).
OGGM is a flowline model that simplifies the glacier geometry by representing it as lines that depict the central flow path.The flowlines are defined following the approach described by Kienholz et al., 2014.The model employs the isothermal shallow ice approximation, assuming that the ice thickness is small compared to its lateral extent, meaning that x-derivatives of stress and velocity are small compared with the z-derivatives (Paterson, 2000).It is important to note that the shallow ice approximation is primarily intended for large ice shelves and may not fully capture the complexities of small mountain glaciers like Monte Perdido.For this glacier, it would be more appropriate to use a model that solves the complete Stokes system, accounting for the three-dimensional nature of ice flow.However, for the sake of simplicity and computational efficiency, we opted to use the OGGM model for this study.Despite this limitation, the validation process supports its use for our specific objectives.
In addition, we used QGIS and CloudCompare software, both open-source platforms, to acquire the outlines of the glacier derived from TLS and compare glacier surface differences between the TLS and OGGM model.

Glacier observation dataset
The surface of the Monte Perdido glacier was derived from terrestrial laser scanning (TLS, RIEGL LPM-321), following the methodology described by López- Moreno et al., (2016).This device generates a 3D point cloud by measuring the distance to thousands of points of the target area with LiDAR technology (Revuelto et al., 2014).TLS observations from 2011 to 2022, allowed us to diagnose the current state of the Monte Perdido glacier and understand its recent evolution (López-Moreno et al., 2019).By analysing the TLS data acquired over this period, we were able to track yearly changes in the glacier's surface elevation, thickness, and extent, providing crucial information about its dynamic behaviour.
The high-resolution topography of the glacier's surrounding area was obtained from the Centro Nacional de Información Geográfica (CNIG) (CNIG, 2023) (Fig. 2).Specifically, we utilised the Digital Elevation Model (DEM) with a 5 m grid resolution.This data allowed us to accurately represent the terrain and its influence on the glacier's behaviour.The glacier outlines were derived from the 2011 TLS (first year with observations), combining this information with topographical variables (Revuelto et al., 2022).Additionally, GPR measurements were obtained in 2016 to capture the ice thickness of the Monte Perdido glacier along several observation transects with an uncertainty of 5 m (López- Moreno et al., 2019).
We specifically focused on analysing the lower Monte Perdido glacier due to the higher availability of TLS and GPR (ground-penetrating radar) data.

Climate data
We used long-term  monthly mean temperature and precipitation data obtained from a weather station located at the Góriz refuge (42.66335°N, 0.01501°E).This meteorological data, managed by the Spanish Meteorological Service (AEMET), was collected approximately 2.5 km from the glacier and at an elevation of 2195 m a.s.l.To adapt the data to the glacier region, we applied a lapse rate of 6.5 °C/km and a precipitation correction factor, which will be further explained in detail.
To ensure a continuous climate dataset, we addressed missing data from the weather station by utilising ERA5 reanalysis data (Hersbach et al., 2023).Prior to filling the gaps, we evaluated the Góriz and ERA5 data during overlapping periods and applied the appropriate multiplication factor to precipitation and temperature.This approach was necessary to maintain data continuity for the OGGM.Furthermore, we have also used the GSWP3-W5E5 data set (Lange and Büchner, 2020), which is the default OGGM climate data.GSWP3-W5E5 dataset is a merge between the GSWP3 (Global Soil Wetness Projected phase 3) dataset (Dirmeyer et al., 2006;Kim et al., 2017) and the W5E5 (biasadjusted ERA5 reanalysis) dataset (Lange, 2019;Cucchi et al., 2020) at 0.5°x 0.5° spatial resolution.
For future projections, we selected a list of 10 global climate models from the Coupled Model Intercomparison Project Phase 5 (CMIP5) ensemble (Taylor et al., 2012) (Table 1).By considering multiple GCMs, we aimed to capture a range of potential future climate scenarios and assess their impact on the glacier.
For CMIP5, four Representative Concentration Pathways (RCPs) have been formulated that provide insights into the expected levels of radiative forcing in the year 2100 when compared to preindustrial conditions (Taylor et al., 2012).These pathways serve as estimations for the impact of greenhouse gas concentrations on the Earth's energy balance.Radiative forcing represents the net change in the energy balance of the Earth system, determined at the top of the atmosphere or the tropopause, due to natural or human-induced perturbations (Myhre et al., 2013).For our analysis, we focused on the RCP2.6 and RCP8.5 scenarios, as they represent the two extremes within the RCP framework.RCP8.5 stands as the "high" scenario, projecting a continuous rise in radiative forcing throughout the twenty-first century until it reaches approximately 8.5 W/m 2 by the end of the century.Conversely, the RCP2.6 scenario assumes strong mitigation efforts with a radiative forcing of 2.6 W/m 2 by 2100 (Taylor et al., 2012).

Model calibration
To determine the volume evolution of our glacier, it is important to know how melt in a given period relates to the climate in the same period.This relationship is established by analysing a period in which we have available both climate data and thickness change, i.e., from 2011 to 2022.Once this calibration is done, we can predict glacier evolution applying future climate forcing to our glacier, assuming that the glacier response to climate forcing remains constant in the future.
For this calibration process, several steps using OGGM were involved.Firstly, we set up the geographical input data for the glacier, such as outlines and local topography.Then the climate data was processed from a user-defined climate file, and later the glacier flowlines were determined (Maussion et al., 2019).Afterwards, we proceeded with the mass balance calibration process.For this, we employed OGGM's standard mass-balance (MB) model, which utilises a temperature index approach (Maussion et al., 2019;Vlug, 2021).The basic assumption of these models is that the melt is proportional to the positive temperature in a certain period of time (Braithwaite and Zhang, 2000;Hock, 2003).This calibration determines specific glacier simulation parameters: the temperature bias, the precipitation factor and the degree-day factor (Maussion et al., 2019;Schuster et al., 2023).
The monthly temperature index model can be calibrated on any mass balance product.The default is the geodetic MB data from Hugonnet et al., 2021, which consists of comparing the glacier surface, obtained from satellite elevation datasets, over two dates (Belart, 2018).This global geodetic glacier dataset provides a mean specific glacier MB estimate between 2000 and 2019 for almost every glacier on Earth (more than 200,000) (Schuster et al., 2023).However, these geodetic estimates do not capture interannual variations and its spatial resolution is moderate when compared to TLS data.Alternatively, in situ mass balance measurements can be employed to capture the year-to-year changes.
The monthly mass balance   at elevation  is computed as follows: Where monthly solid precipitation    is multiplied by the precipitation correction factor  .As there is no precipitation lapse rate in the model,   can be seen as a global correction factor for orographic precipitation, avalanches, and wind-blown snow (Vlug, 2021).The precipitation is assumed as liquid above 2°C, solid below 0°C, and the fraction of solid precipitation is linearly interpolated between these two boundary values.  is the monthly mean air temperature at 2 m and   is the monthly mean air temperature above which ice melt is assumed to occur (-1°C per default according to OGGM standard due to ice pressure).The temperature lapse rate is set by default to 6.5 °C/km.The parameter   is the degree-day factor indicating the temperature sensitivity of the glacier (van der Laan et al., 2022;Schuster et al., 2023).

PAPER ACCEPTED. PREPRINT VERSION
We conducted three simulations (each with a specific calibration of three parameters: precipitation factor, temperature bias, and degree-day factor) to analyse the behaviour of glaciers under different configurations (Fig. 3, Table 2):

PAPER ACCEPTED. PREPRINT VERSION
-Simulation 1: This simulation involved calibrating the mass balance using in situ data obtained from TLS for the years 2011 to 2022.The precipitation factor was estimated by comparing the total water equivalent of snow during the accumulation period (from October to April) at the glacier location (López- Moreno et al., 2019) with the precipitation data recorded at the Góriz weather station for the same period.Analysis of specific years (2013-14, 2014-15, and 2016-17, which are the only periods with glacier surface observations during the accumulation period) revealed a maximum mean snow accumulation in late April of 3.25 m, with an average snow density of 454 kg/m 3 (López- Moreno et al., 2019).This indicated a total water equivalent of 1475.5 mm for the accumulation period.Considering that the mean precipitation observed at the Góriz weather station during the same period was 1089.6 mm, a precipitation factor of 1.3 was estimated and used for the mass balance calibration.Then, the model adjusted the degreeday factor and the temperature bias to minimise the difference between the model outputs and observed data, ensuring a better fit between the simulated and the actual mass balance.
-Simulation 2: As in the first simulation, this one involved mass balance calibration using TLS data.However, instead of using weather station climate data, we used the default climate data (GSWP3-W5E5) provided by the OGGM framework.The three parameters (precipitation factor, temperature bias, and degree-day factor) were adjusted accordingly.
-Simulation 3: In this simulation, the mass balance was calibrated with the default average geodetic observations from January 2000 to January 2019 of Hugonnet et al., 2021, and the default climate data (GSWP3-W5E5) provided by OGGM.
All the calibration alternatives have at least the three free parameters mentioned above (Schuster et al., 2023).Without these parameters, the observed glacier MB often cannot be reproduced by the model.
Once the mass balance calibration is performed, a standard geometry evolution model, which is a depth-integrated flowline model, is responsible to compute the change in glacier geometry.Before running this simulation, stable glacier conditions are required at the beginning of the study period.To ensure this, we performed a spin-up process, where the geometry and evolution of the glaciers were initialised from a given year.We selected a fixed geometry spin-up year of 2000, approximately 10 years before the date of the outline (2011).This year (2011) represents the point at which the glacier is expected to reach equilibrium (Maussion et al., 2019).
Using the flowline model, an estimate of the ice flux along each glacier grid point cross-section is computed by making assumptions about the shape of the cross-section (parabolic, rectangular or trapezoid) and relying on mass-conservation consideration (Maussion et al., 2019).Using the physics of ice flow and the shallow ice approximation, the model then computes the thickness of the glacier along the flowlines and the total volume of the glacier.
After performing the historical climate run, we used the ten GCMs described in Table 1 to project future volume of the Monte Perdido glacier.We employed an ensemble of GCMs to represent the temperature and precipitation variability in climate projections.To downscale global climate data to a regional level, we obtained precipitation and temperature data for each GCM from the OGGM server hosted by the University of Bremen and we applied the precipitation and temperature biases, to the baseline local climatology, which was defined using the Góriz weather station climate dataset.The model then uses the interpolated climate data for each GCM to calculate glacier mass balance.
Using this climate data, we ran the simulation for each GCM and scenario from 2020 to 2100.To initiate this task, we inputted the GCM climate data and utilised the spun up geometry and mass balance conditions from the historical run as initial conditions within the model.

PAPER ACCEPTED. PREPRINT VERSION
Finally, we compiled the 20 simulations generated from the ensemble of GCMs and merged them into two datasets, one for each RCP scenario.This allowed us to calculate the median values and plot the evolution of glacier volume (Fig. 7).

Evaluation of climate data
We first evaluate the GSWP3-W5E5 dataset from the OGGM repository against the data measured at the Góriz weather stations.
Figure 4a represents the historical annual temperature data measured at the Góriz weather station, along with the GSWP3-W5E5 dataset.The data reveals long-term temperature variations in the region, and a notable temperature rise on both datasets.The 30-year rolling average highlights this increase smoothing out short-term variations.The slope of 0.04°C/year and the p-value < 0.05 on both datasets, confirm a statistically significant positive trend in temperature.

PAPER ACCEPTED. PREPRINT VERSION
In addition, Figure 4b shows the historical total annual precipitation data obtained from the weather station and the GSWP3-W5E5 dataset.The data provides the interannual variability of precipitation over time.Unlike the temperature, there is no significant change in precipitation over time (p-value > 0.05).Furthermore, we observed a mean temperature difference of 0.4°C between the Góriz weather station, situated at 2195 m a.s.l., and the GSWP3-W5E5 dataset, located at 1756 m a.s.l.Additionally, the Góriz weather station records a greater total annual precipitation compared to the GSWP3-W5E5 dataset, with a difference of 590 mm.

Comparison of modelled and TLS-derived volume differences
The comparison of TLS volume evolution (surface decrease multiplied by glacier extent) and modelled volume evolution provides insights into the accuracy of the model in replicating the observed glacier volume changes (Fig. 5).The root mean square error (RMSE), correlation coefficient, and pvalue were calculated to evaluate the model performance (Table 3).The results indicate that both simulations with in situ mass balance calibration, namely simulation 1 and 2, exhibit a higher correlation with the TLS-derived volume compared to simulation 3.This suggests that the calibration process with TLS data improves the agreement between the model and the observed data.Furthermore, simulation 1 with Góriz climate data shows slightly improved performance compared to the one with GSWP3-W5E5 climate data, as evidenced by a lower RMSE value and higher correlation coefficient (Table 3).

PAPER ACCEPTED. PREPRINT VERSION
The incorporation of TLS data into the calibration process helps to account for the interannual variation of the glacier mass balance.This, in turn, enhances the accuracy of the model's predictions.Furthermore, Góriz climate data might provide a better representation of the local climate conditions and their influence on the glacier, resulting in a more accurate estimation of mass balance parameters.
It is important to note that simulation 3, despite exhibiting a slightly lower correlation and higher RMSE, still demonstrates a reasonable agreement with the TLS derived volume.This suggests that the calibration process, even with GSWP3-W5E5 climate data and geodetic mass balance calibration, can provide valuable insights into the glacier's behaviour.However, the differences in performance between simulation 3 and the other simulations highlight the importance of selecting appropriate climate data and calibration methods to improve the accuracy of glacier volume projections.

Comparison of modelled and GPR thickness
The GPR measurements, taken in 2016, offer a direct assessment of the glacier's ice thickness at specific locations on the Monte Perdido glacier (López- Moreno et al., 2019).The comparison of modelled and GPR thickness provides additional insights into the accuracy of the model's representation of the spatial distribution of ice (Fig. 6).It has to be noted that the GPR measurements were influenced by the presence of water, causing a very low signal to noise ratio leading to ± 5 m of uncertainty in estimations of ice thickness, which should be taken into account when interpreting the results.
The mean difference between the modelled and GPR thickness is 6.4 m, with a maximum difference of 32.7 m, indicating a certain level of variability and uncertainty in the model's ability to capture the ice thickness distribution.Furthermore, the average discrepancy between the two datasets is provided by the RMSE of 8.7 m.While there is some level of agreement between the model and the GPR data, there are still considerable differences between them.
However, it is important to note that GPR measurements provide localised information and may not fully represent the entire glacier's ice thickness distribution.Additionally, accurately modelling ice thickness is challenging without precise knowledge of the topography below the glacier.

PAPER ACCEPTED. PREPRINT VERSION
Moreover, the weather station used for collecting meteorological data for the Monte Perdido glacier is situated about 2.5 km away from the glacier and on the south face, while the glacier itself is on the north face.This spatial difference introduces uncertainty in representing local climate conditions, despite applying temperature and precipitation correction factors.The use of a fixed lapse rate of 6.5°C/km to adjust weather station data to the glacier region simplifies temperature variations with elevation.Additionally, the estimation of the precipitation factor based on comparing data from the glacier location and a weather station may overlook variations in snowfall patterns, snow density, or liquid precipitation during the accumulation period.

Future volume projections
Across both RCP scenarios, the projected total ice volume for the Monte Perdido glacier shows a consistent decrease from 2020 to 2100. Figure 7 shows a faster decline in volume between 2020 and 2060, followed by a deceleration in the rate of decrease.

PAPER ACCEPTED. PREPRINT VERSION
The total median volume exhibits minimal variation between the RCP scenarios.For the RCP8.5 scenario, the ice volume experiences a significant decrease of 95.2% by the year 2100, 88.6% by 2060, and 66.0% by 2040.Similarly, under the RCP2.6 scenario, there is a decrease of 91.6% by 2100, 82.8% by 2060, and 59.6% by 2040.
The observed decreasing trend in the volume of the Monte Perdido glacier is not unexpected; many studies have documented similar decreases in glacier volume worldwide (Ma et al., 2010;Zekollari et al., 2019;Khadka et al., 2020), including in the Pyrenees (Chueca Cía et al., 2005;López-Moreno et al., 2016;Campos et al., 2021;Vidaller et al., 2021).Given that Monte Perdido is one of the largest glaciers of the Pyrenees, situated at higher altitudes and facing north, it suggests that other glaciers might experience even more pronounced losses.

Conclusions
We employed the Open Global Glacier Model to simulate and analyse the Monte Perdido glacier recent and future evolution.We utilised the OGGM to estimate the total ice volume of the glacier and simulate its evolution in response to different climate forcings.
The calibration process of the OGGM involved three simulations: one using in situ mass balance calibration and weather station climate data, another with in situ mass balance calibration using default climate data, and a third with uncalibrated mass balance and default climate data.Through these simulations, we evaluated the performance of the model when replicating the observed volume changes of the glacier.The simulations with in situ mass balance calibration exhibited a higher correlation with TLS-derived volume compared to the default MB geodetic calibration, indicating the importance of exploiting in-situ observations on the calibration to improve model accuracy.
Furthermore, we projected the future volume of the Monte Perdido glacier using an ensemble of ten GCMs under the RCP2.6 and RCP8.5 scenarios.The results showed a consistent decrease in total ice volume from 2020 to 2100, with a faster decline between 2020 and 2060 followed by a deceleration in the rate of decrease.The projected volume reductions were substantial, ranging from 91.6% to 95.2% by the year 2100, depending on the scenario.These findings align with the global trend of glacier volume decrease and are consistent with previous studies in the Pyrenees region.
In addition, as we look ahead, it's worth considering future directions that could further enhance our understanding of glacier behaviour and its interaction with climate.While our study primarily focused on the Monte Perdido glacier's response to climate forcing, we acknowledge the need for future investigations into the potential consequences of climate change on avalanche triggering and its subsequent impact on glacier evolution.Moreover, new approaches like ODINN.jl (Bolibar et al., 2023) and MuSA (Alonso-González et al., 2022) demonstrate promising paths to enhance glacier modelling by incorporating advanced data assimilation techniques.

Figure 1 :
Figure 1: (a) Location and extent of the Monte Perdido glacier in 2022 (coordinates in extended UTM zone 31 T).(b) View of Monte Perdido glacier on October 5, 2022.

Figure 3 :
Figure 3: Comparison of TLS-derived mass balance (and its standard deviation (SD) and (a) modelled mass balance with in situ MB calibration and Góriz weather station climate data, (b) modelled mass balance with in situ MB calibration and GSWP3-W5E5 climate data, and (c) modelled mass balance with geodetic MB calibration.

Figure 4 :
Figure 4: (a) Annual mean temperature and (b) total annual precipitation obtained from Góriz weather station (grey) and GSWP3-W5E5 dataset (green).The 30-year rolling average is depicted in dashed lines.

Figure 5 :
Figure 5: Monte Perdido glacier volume (and SD) since 2011 derived from TLS data compared with the three simulations.

Figure 6 :
Figure 6: Difference between the modelled ice thickness and the GPR measurements taken in 2016 at specific locations of Monte Perdido glacier.WGS84 coordinate system.

Figure 7 :
Figure 7: Multi-GCM ice volume for RCP2.6 and RCP8.5 scenarios from 2020 to 2100 with the historical simulation from 2000 to 2020 in a dashed line and the calibrated period in a solid black line.Volumes for each GCM run of RCP2.6 are plotted in blue and for RCP8.5 in red, with multi-GCM medians represented in thicker lines.Plot (a) displays the full projection, while (b) zooms the 2020 -2080 period.

Table 1 .
CMIP5 models used in this study.

Table 2
. Data sources used for mass balance calibration and climate data in different simulations

Table 3 .
Comparison of TLS-Derived Volume with the three simulations from 2011 to 2020