Seasonal water storage and release dynamics of bofedal wetlands in the Central Andes

Tropical high‐Andean wetlands, locally known as ‘bofedales’, are key ecosystems sustaining biodiversity, carbon sequestration, water provision and livestock farming. Bofedales' contribution to dry season baseflows and sustaining water quality is crucial for downstream water security. The sensitivity of bofedales to climatic and anthropogenic disturbances is therefore of growing concern for watershed management. This study aims to understand seasonal water storage and release characteristics of bofedales by combining remote sensing analysis and ground‐based monitoring for the wet and dry seasons of late 2019 to early 2021, using the glacierised Vilcanota‐Urubamba basin (Southern Peru) as a case study. A network of five ultrasound loggers was installed to obtain discharge and water table data from bofedal sites across two headwater catchments. The seasonal extent of bofedales was mapped by applying a supervised machine learning model using Random Forest on imagery from Sentinel‐2 and NASADEM. We identified high seasonal variability in bofedal area with a total of 3.5% and 10.6% of each catchment area, respectively, at the end of the dry season (2020), which increased to 15.1% and 16.9%, respectively, at the end of the following wet season (2021). The hydrological observations and bofedal maps were combined into a hydrological conceptual model to estimate the storage and release characteristics of the bofedales, and their contribution to runoff at the catchment scale. Estimated lag times between 1 and 32 days indicate a prolonged bofedal flow contribution throughout the dry season (about 74% of total flow). Thus, our results suggest that bofedales provide substantial contribution to dry season baseflow, water flow regulation and storage. These findings highlight the importance of including bofedales in local water management strategies and adaptation interventions including nature‐based solutions that seek to support long‐term water security in seasonally dry and rapidly changing Andean catchments.


| INTRODUCTION
In the Southern Tropical Andes, high-altitude wetlands are a characteristic feature of the mountain hydrology. They are locally known as 'bofedales' or 'oconales' (Chimner et al., 2019;Gandarillas et al., 2016;Garcia & Otto, 2015) and cover most parts of Central and Southern Peru, Bolivia, Northern Chile and Argentina at an altitudinal range of about 3500-5200 m asl (García & Beck, 2006;Meneses et al., 2019;Otto et al., 2011;Squeo et al., 2006). These ecosystems are characterized for being remarkably exposed to extreme conditions, such as strong winds, diurnal thermal contrasts, pronounced dry seasons, high solar radiation, and low oxygen concentrations (Gandarillas et al., 2016;Squeo et al., 2006). Bofedales are known to be high primary production systems of seasonally dry and moist habitats (Buitr on-Aliaga & Fernández-Callisaya, 2012), typically located on flat or slightly inclined (from 0 to 15 ) fluvial and glacial valleys, alluvial plains and fluvial terraces (Garcia Dulanto, 2018;Mango-Mamani, 2017). They develop over colluvial deposits or plane surfaces, such as lacustrine plains, or filled shallow pits, such as fluvioglacial kettles or bed-rock potholes (Squeo et al., 2006) and are inhabited by a diverse set of low-growing hydrophytic plants that often accumulate peat and are occasionally accompanied by herbaceous vegetation between 0.1 and 0.5 m tall (MINAM, 2019).
Generating a well-founded understanding of high-Andean hydrology, including in wetland dominated regions, has been difficult due to a strong lack of resources and hydrometeorological monitoring. Findings are fragmented as catchment features (i.e., glaciers, wetlands and groundwater) tend to be researched in isolation (Navarro et al., 2023).
Few studies have aimed to characterize the role of bofedales in catchment hydrology (Cooper et al., 2019;Polk et al., 2017;Valois et al., 2020;Valois et al., 2021), however, findings are widespread and somewhat inconsistent. Furthermore, no attempts have been made to quantify their relative streamflow contribution or seasonal buffering capacity. This evidence would be paramount in building an understanding of how bofedales will respond to changes in climate and human intervention and, therefore, how future downstream water supply might be affected. This lack of hydrological characterization combined with the vulnerability of bofedales to human disturbance and environmental change, and their importance for local and regional water supply, poses severe challenges to their sustainable management (e.g., Drenkhan et al., 2023). In this context, using the headwaters of the Vilcanota-Urubamba basin (Southern Peru) as a case study, we provide a framework to quantify the seasonal water storage and release dynamics of bofedales using minimal monitoring and widely available data. This study aims to: (1) understand the spatio-temporal dynamics of bofedales extent by mapping them in respective wet and dry seasons, and (2) estimate the streamflow contribution of bofedales relative to other catchment features, using a novel combination of in situ measurements, analysis of remotely sensed data, and hydrological modelling.

| Study region
Our study catchments are located within the Vilcanota-Urubamba basin (VUB) in Southern Peru in the transition area of the wet Puna grassland and steppe ecoregion (Squeo et al., 2006) in the south and the tropical rainforests in the north of the basin (Figure 1). The VUB covers a surface area of 11 052 km 2 and an elevation range from about 1100 to 6300 m asl (Figure 1a,b). It is a glacierised basin with total glacier area of 142 km 2 (Drenkhan et al., 2018). Our study focuses on the Sibina Sallma (41.6 km 2 ) and Halairipampa (54.1 km 2 ) catchments ( Figure 1c) in the high wet Puna region in the southeastern part of the VUB. Halairipampa includes a large glacier coverage (about 8.9 km 2 ), whereas glacier extent in Sibina Sallma is very limited (<0.1 km 2 ).
Two weather stations are in the vicinity of the studied catchments: Sibinacocha Lake and Quelccaya Glacier (4880 and 5650 m asl, respectively; Figure 1c). Sibinacocha indicates 776 mm/year precipitation on average (2017-2021) from which 66% (511 mm) occurred during the wet season (December-March) and 4% (29 mm) during the dry season (June-August). Average annual precipitation at Quelccaya was 1199 mm/year (2016-2020) from which 60% (718 mm) occurred during the wet season and 7% (79 mm) during the dry season. Annual average surface temperature at Sibinacocha varied between 2.6 C in austral wet summer (December-March) and 1.0 C in austral dry winter (June-August).
Exposed bedrock on the Sibina Sallma valley slopes is dominated by the Devonian metapelitic Cabanillas group which is locally overlain by the younger volcano-sedimentary Mitu group to the south. The Halairipampa catchment bedrock is dominated by the Carboniferous Ambo (sandstone) Group. The bofedales in both catchments are almost exclusively constrained to the valley floors where they overlie extensive, highly permeable glacio-fluvial and alluvial deposits. They typically comprise a high-permeability organic peat layer underlain by clay-rich lacustrine deposit. This low permeability horizon likelypromotes the formation of the bofedales by inhibiting surface drainage to the buried unconsolidated deposits beneath (Glas et al., 2018).
Total bofedal cover for Peru has been estimated by the Peruvian Ministry of Environment (MINAM) within the national ecosystem inventory with a baseline for 2016. According to this assessment, the Sibina Sallma catchment covers a bofedal area of 3.8 km 2 (9.1%) while Halairipampa includes 5.8 km 2 (10.7%), bofedal area over the VUB is 267.4 km 2 (2.4%) and 5481.7 km 2 (0.4%) for Peru (MINAM, 2019). This inventory, partially based on remote sensing data and only considering larger bofedal areas (minimum = 20 ha), has important limitations as total bofedal area may highly depend on seasonal variability, and the considered thresholds and model accuracy (UNEP, 2022).

| Seasonal bofedales extent mapping
Bofedales is a local term used across the Southern Tropical Andes that does not yet have a formal definition. For example, in some literature, bofedales are considered peatlands (e.g., Valois et al., 2020;Valois et al., 2021). However, we have adapted a version of the definition given in Maldonado Fonkén (2014) to areas of land that are saturated for a sufficient period during the year to sustain various types of highaltitude wetland plant communities. The underlying soil may or may not have the presence of peat, they can be seasonally or permanently saturated and natural or man-made. Therefore, our classification model aims to capture all high-altitude wetlands within the study region, including but not limited to peatlands. The model was used to map bofedales extent in the dry and the wet season. However, we do not assume that actual bofedal area changes seasonally but, rather, that the resulting maps give an indication to the seasonal differences in bofedal water storage.
Bofedales extent for the VUB was mapped using multi-spectral Sentinel 2 observations. This dataset is particularly suited for vegetation analysis because it includes a near-infrared band at 10 m as well as three red-edge bands and two short-wave infrared bands at 20 m spatial resolution (Gatti & Bertolini, 2015). We acquired Sentinel 2 observations for the end of the dry (01/08/2020-30/11/2020) and wet (01/04/2021-14/06/2021) seasons. Images were filtered by pixel-wise cloud cover probability >10%, then mosaicked into cloudfree images for each season using median band values (Shepherd et al., 2020).
Next, we applied supervised classification using the Random Forest (RF) machine learning algorithm. This method has been proven to provide uniform and high accuracy results when applied over large areas (Mutanga & Kumar, 2019) and has been successfully applied to map tropical high-Andean wetlands (cf. Chimner et al., 2019;Hribljan et al., 2017;Jara et al., 2019). RF uses bootstrap aggregation (bagging) to generate an ensemble of 'trees' in which each classifier is trained on a random subspace of features (Belgiu & Dr aguţ, 2016). Once trained, the final prediction is obtained by a majority voting system of each model from the ensemble (Ham et al., 2005;Mutanga & Kumar, 2019).
A major benefit of the RF method is the ability to incorporate covariables to improve the classification (Belgiu & Dr aguţ, 2016;Ham et al., 2005). Here, we used a total of 20 vegetation and 7 topographic indices (Appendix A, Tables A1, A2) that potentially enhance the spectral differentiation of vegetation from satellite imagery to identify potential bofedales extent (Jara et al., 2019;Mahdavi et al., 2018).
Topographic indices were derived from the NASADEM dataset, which has a 30 m spatial resolution (Buckley et al., 2020).
In absence of ground-validated data on land cover types within the study area, a training point sampling strategy was developed using two Regions of Interest of clearly identifiable bofedal and non-bofedal areas. Therefore, iteratively determined thresholds were defined based on a scoping literature review using the three most common vegetation and topographic indices for wetland classification for both the dry and wet season (Appendix A, Tables A3-A5).
A stratified random sampling method was applied to the resulting two Regions of Interest to reduce the accuracy estimation errors of non-abundant features (i.e., bofedal pixels) and to increase the number of samples to improve the accuracy of the model without biasing the performance estimators (Olofsson et al., 2014). The sampling size necessary to achieve a target standard error was calculated following Cochran (1977). The RF model was then refined to reduce overfitting (Appendix A, Random Forest (RF) model refinement).
To evaluate model performance, a binary confusion matrix was generated based on the Matheus Correlation Coefficient (Chicco et al., 2021), sensitivity and accuracy (Starovoitov & Golub, 2020) performance assessments. Lastly, three postprocessing steps were applied to reduce residual errors. The first two were adopted in Baum & Zanetti (2015), Kumar (2020) and Mohammadi et al. (2020): (i) eliminating bofedal pixels if less than half of their eight neighbour pixels were considered bofedales and (ii) eliminating bofedal areas less than 3 Â 3 pixels. Pixel resolution is attributed to the 30 m resolution NASADEM input data. The outcome resolution encompasses a minimum of 3 Â 3 pixels (0.8 ha) necessary in determining a category to map. This resolution has been shown to be effective to identify vegetation (Li et al., 2021;Shamsoddini & Raval, 2018) and, specifically, bofedales (Garcia Dulanto, 2018;García & Lleellish, 2012). The final postprocessing step, (iii) a water mask based on the Normalized Difference Water Index was applied to avoid confusion of bofedales with lakes and other open water bodies.

| Hydrometeorological data collection and processing
Groundwater observation boreholes of 0.15 m diameter were installed at depths between 1.1 and 2.5 m in four representative bofedal locations across the two catchments (US1-4, Figure 1; Table 1).
For boreholes US2 and US3, total borehole depth was limited upon reaching hard clay lacustrine material where manual core extraction was not feasible. Afterwards, the borehole was widened to allow the insertion of perforated pipes. The pipes were capped with a PVC lid onto which a downward-facing ultrasound water level logger was installed. An additional ultrasound logger (US5) was fitted at the outlet of the Halairipampa catchment to measure river levels. At this location, previous water level data from a pressure transducer (P1) are also available. The records between US5 and P1 showed good agreement for the overlapping period and, therefore, they were merged into a single time series of water level. US1-5 recorded data at a 5-min interval (P1 at a 30-min interval).
The measurements were datum-adjusted and smoothed using a moving average with a window of 30 min to reduce noise. These were averaged to daily frequency to minimize the impact of the diurnal temperature cycle in these data. This is caused by temperature variations in the air column of the pipe, which are exceptionally high in tropical mountain regions and affect ultrasound wave transmission.
Three manual flow estimations were made between October 2016 and December 2018 at the outlet of the Halairipampa catchment using float gauging. A power law was fitted to the data points to formulate a rating curve. The resulting discharge hydrograph was separated into stormflow and baseflow following Gustard et al. (1992). We assume that all stormflow is produced by surface runoff due to the velocity of the response. deep groundwater contributions to hillslope or wetland recharge and streamflow are often negligible in high mountain conditions. This may not be the case in our catchment because of the presence of proglacial moraine deposits which act as groundwater stores (Somers & McKenzie, 2020). However, it is a conservative assumption because it will lead to an underestimation of the catchment's hydrological buffering capacity.
Conceptual models to estimate each of these flows from our observations are outlined in the next subsections. Because of a lack of reliable radiation data, we estimated evapotranspiration by closing the catchment water balance using precipitation and discharge data.
Lastly, it is worth pointing out that isotope analysis by Cooper et al. (2019) found that streamflow and glacier melt contributions to hillslope or wetland recharge in this region were negligible. Therefore, both are modelled as parallel and independent flow pathways. We apply our model to the Halairipampa catchment, since only discharge data for this catchment were available. Resolving the catchment water balance for the 16-month study period using the Sibinacocha rainfall data, we found a runoff ratio of 1.00 for Halairipampa, with total rainfall and runoff equal to 1371 mm. This is unrealistically T A B L E 1 Characteristics of the hydrological loggers used for this study in the two sub-catchments Sibina Sallma and Halairipampa. high and indicates errors in the water balance. We attribute this error to the location of the Sibinacocha meteorological station, which is southwest of the Halairipampa catchment. A strong northeastsouthwest gradient of decreasing precipitation is present in the region because of moisture spill-over from the Amazon basin. The Quelccaya rainfall data, which gave a runoff ratio of 0.897 for the studied period, were used as the driving input to the model. The runoff ratio is still unrealistically high; however, our interests are in the relative streamflow contributions of end members and their temporal dynamics, rather than predicting flow magnitudes. As such, we did not apply a correction factor.

| Surface runoff (Q S )
Studies suggest that surface runoff production in high-Andean catchments is predominantly caused by saturation excess rather than infiltration excess (Buytaert & Beven, 2011;Lazo et al., 2019;Mosquera et al., 2015;Mosquera et al., 2016). Therefore, we posited a direct linear relation between bofedales extent and surface runoff production in the catchment. However, a linear regression between the monthly surface runoff percentage of rainfall and the monthly average areal bofedales extent gave a weak negative correlation ( p-value = 0.823).
Therefore, a constant factor (total surface runoff over total precipitation) was used to simulate daily surface runoff as a fraction of precipitation.

| Shallow subsurface runoff (Q W+H )
Previous studies have shown that the subsurface runoff of high-Andean wetland catchments behaves as a set of parallel linear reservoirs with different residence times (Buytaert et al., 2004;Buytaert & Beven, 2011). Here we combine this linear reservoir behaviour with the observational data from the water level sensors to build a conceptual model of the shallow subsurface runoff. We model Q W and Q H jointly as Q W+H because hillslopes often feed into wetlands and our experimental approach does not allow separation of the fluxes. A linear reservoir is characterized by the equation: where Q is the outflow of the wetlands, k is the linear reservoir drainage rate coefficient (rate constant) and S is the total storage capacity of the wetlands and hillslopes. If we assume that wetland storage is much larger than hillslope storage and that there is a linear relation between the wetland extent and the storage capacity, then S can be represented as: where A w is the maximum wetland extent during the wet season, A d is the minimum wetland extent during the dry season, h is the water level of the wetland above the seasonal minimum, and H is the maximum level difference between A d and A w . Combining both formulas yields: where u = unit hydrograph ordinate, v = volume of direct runoff, t = time, k = parameter with dimension of time, n = dimensionless parameter and Γ = gamma function.
First, the lag time between precipitation input and the unit hydrograph peak was determined by maximizing the cross-correlation between rainfall and subsurface water levels (Cooper et al., 2019). Lag times for bofedales in the Sibina Sallma catchment were calculated using data from the Sibinacocha weather station, based on their proximity, whereas lag times for Halairipampa were calculated using data from the Quelccaya weather station (Figure 1c). During the dry season, levels at US3 fall below the depth of the well, therefore, this section of the data was removed from the cross-correlation analysis.
A smoothing spline was applied to the US2 and US4 cross-correlation plots to reduce noise.
For the rising limb, a linear increase to the peak was assumed. For the falling limb, respective sections of the observed recession curve during the dry season, least disturbed by precipitation, were used (Appendix C, Table C1). These were extended linearly to a point of zero flow and the results were normalized to unit volume. The Nash (1957) distribution was fitted to these unit hydrographs using the Nelder-Mead algorithm (Nelder & Mead, 1965) to minimize the rootmean-square error (RMSE) with the Nash (1957) parameters k and n, and the point of zero flow as free parameters. Since the US2 and US4 unit hydrographs showed high resemblance, a combined unit hydrograph was fitted to the data from both locations. This resulted in 3 unit hydrographs. As we selected representative points in the landscape, we assume that these jointly represent the bofedales behaviour at the catchment scale.
Subsequently, to estimate the relative contribution of each of those hydrographs to the catchment flow, we introduce partitioning constants to weigh each hydrograph, such that: where P is the precipitation time series, U i are the unit hydrographs, * is the convolution operator, and w i are the partitioning constants, which are greater than 0 and sum to 1.
The values for these partitioning constants were estimated by minimizing the RMSE between the observed subsurface flow (baseflow minus glacial melt) and modelled subsurface flow, again using the Nelder-Mead algorithm (Nelder & Mead, 1965). There is a distinct reduction in the Quelccaya 2020-2021 wet season rainfall which is not reflected in the Halairipampa hydrograph. Therefore, we fitted the modelled subsurface flow to the observed 2019-2020 wet season and following dry season subsurface flow (until 11/08/2020).
Lastly, we assume that the outflow of the bofedales contributes directly to river discharge without any further delay, for example by passing through groundwater storage. This is a conservative assumption because it will likely lead to an underestimation of the bofedales' buffering capacity.

| Bofedales extent
The supervised classification using the RF algorithm reveals a total bofedal area of 1.4 km 2 and 6.3 km 2 for the Sibina Sallma and Halairipampa catchments, respectively, during the dry season (2020). During the wet season (2021)  and drying of small lakes has previously been described (Drenkhan et al., 2018). However, the results presented in this study are limited by the short study period focusing on one specific dry and wet season in the periods 2020-2021. More research is needed to understand the seasonal and interannual interplay of local climate (e.g., changes in the precipitation regime), morphology and deposits (e.g., stratification of permeable and impermeable layers), the role of meltwater and other surface and subsurface inputs, and human impacts in the catchments (e.g., bofedal management, irrigation and peat extraction).

| Water level data
All water level data reflect the strong climate seasonality (Figure 4), in line with areal bofedales extent as observed from the satellite imagery. However, some clear differences in bofedal storage behaviour can be observed. This heterogeneity is reflected in lag times obtained from the cross-correlation analysis (Table 2) and is compatible with F I G U R E 4 Water table observations from logger network (see Figure 1c and Table 1 for the exact locations).
the findings of Cooper et al. (2019). This is not surprising because the response times of wetlands will be determined strongly by their position in the landscape, which controls both the upslope contributing area and the velocity of water release. As mentioned in the methods, US3 water levels drop below the depth of the well during the dry season and were removed from further analysis.

| Surface runoff production
We did not find a statistically significant correlation between bofedales extent and surface runoff production ( p-value = 0.823); in fact, the fitted regression has a slightly negative trend. This is counterintuitive but may be caused by the level of hydrological connectivity.
Increasing areas of disconnected wetlands may rather enhance infiltration due to improved hydraulic conductivity thereby reducing surface runoff in the catchment (Buytaert & Beven, 2011). This may act as a confounding variable for the relation, especially if it does not show a clear seasonality. However, the precipitation time series is not representative for the catchment. This is particularly noticeable during the 2020-2021 wet season, which is uncharacteristically low for the Quelccaya weather station and is not reflected in the Halairipampa discharge data ( Figure 6). Furthermore, inconsistencies can be seen in the dry season where the timing of storm events in the rainfall and discharge data do not match. Differences have been attributed to localized climatic heterogeneities, which are prevalent in this region.
Further exploration into the seasonality of surface runoff production would require more representative rainfall data through the installation of rainfall sensors within the catchment.

| Subsurface unit hydrographs
While our approach does not allow separating the hillslope release Q H from the wetland release Q W , the wide range of lag times (Table 2) is indicative for the diverging hydrological response of different bofedales within the catchment. The lag time of 1 day for US1 suggests a wetland that fills quickly, most likely because of overland flow and short hillslopes (Tetzlaff et al., 2014). On the other end of the spectrum US2 and US4 have lag times of around 30 days, which suggest that they are fed by a slower subsurface response, for example from longer hillslopes or deeper subsurface flow (e.g., moraines). To get further insights in the topographic and connectivity controls, we attempted to find a relation between the lag times and topographic indices such as the upslope contributing area and soil topographic index of the bofedal locations, but this was unsuccessful, likely due to the small number of observations and the coarse resolution of the available digital elevation model.
The recession curves also show highly different behaviour ( Figure 5). UH1 and UH3 drain relatively quickly, reaching half of their peak outflow in about 10 days, while UHC only drops to 50% of its peak after more than a month. This suggests that the mechanisms controlling bofedal outflow are also highly heterogeneous and may be related to processes such as the morphology of the bofedales drainage system, as well as the hydraulic conductivity, preferential flow paths and thickness of the soil mounds that separate the bofedales from the streams. Nevertheless, the unit hydrographs provide insight in the temporal dynamics of bofedales' contributions to streamflow, which can remain substantial for several months after precipitation events.

| Combined hydrological model
Calibrating the partitioning constants (Table 2)   F I G U R E 5 Modelled subsurface unit hydrographs fitted to observed lag times and recession curves (see Table 2 for lag times and fitted Nash (1957) (Blöschl & Sivapalan, 1995). Adaptation strategies which prioritize natural infrastructure, such as wetland irrigation (e.g., FAO, 2022a) and restoration (e.g., FAO, 2022b) to buffer seasonally decreasing water quantity need to consider a solid evidence baseline to identify costeffective and hydrologically most beneficial options . For improved relevance in regional, national and international policy, it is therefore essential to develop methods that can quantify the hydrological benefits of bofedales at larger scales. This will require a better understanding of the different types of bofedales to identify representative sites for increased monitoring and subsequent development of regional hydrological models that capture the streamflow contributions of bofedales and their role within the wider basin response.

| CONCLUSIONS
Bofedales may provide critical water storage and seasonal buffering for Andean communities, which are increasingly important in view of the combined impacts of climate and land use changes and growing water demand in downstream areas . In this study, we show how to estimate the storage and release characteristics of bofedales and to understand their role within the wider catchment response using a combination of remote sensing techniques and in situ hydrological monitoring of bofedal water levels, precipitation and streamflow. Applying our methodology to the Vilcanota-Urubamba basin, we find a strong seasonal variation in bofedales' water storage at both the sub-catchment and wider river basin scale. We applied the linear reservoir concept to bofedal water levels to estimate unit hydrographs for representative bofedales. The large differences in the residence time of these unit hydrographs (1-30 days) suggest a strong heterogeneity in subsurface storage and release behaviour between bofedales, which we attribute to hydrological connectivity and upslope contributing area. Upscaling these hydrographs to the catchment scale by means of calibrated contribution factors allowed estimating the contribution of different end members (glacier melt, surface runoff and combined hillslope and bofedal subsurface outflow). We find that nonglacier subsurface flow is substantial (59% of total streamflow over the entire modelling period) and is sustained for several months after precipitation events (74% of core dry season streamflow). This highlights the importance of considering bofedales in local water management strategies including nature-based solutions. It also illustrates the need to base catchment interventions on a solid evidence baseline that identifies the specific hydrological conditions and benefits of high-altitude wetlands to counteract decreasing water availability. Improved monitoring in representative bofedal sites is, therefore, crucial to enable upscaling our hydrological understanding for relevance in regional, national and international policy.

ACKNOWLEDGEMENTS
This study was developed within the framework of the Newton-

CONFLICT OF INTEREST STATEMENT
The authors declare that there is no conflict of interest.

DATA AVAILABILITY STATEMENT
• Glacier and bofedales extent in Figure 1: data available upon request.
• Sibinacocha Lake weather station data: provided by SENAMHI, • Input data and results from the glacial melt model: available upon request.
• Results from the bofedal classification: available upon request.