Surface elevation and ice thickness data between 2012 and 2020 at 2 the ablation area of Artesonraju Glacier , Cordillera Blanca , Perú

We present a representative set of data of interpreted ice thickness and ice surface elevation at the ablation area of 12 the Artesonraju glacier between 2012 and 2020. The ice thickness was obtained by means of Ground Penetrating Radar 13 (GPR), while the surface elevation was by means of automated total stations and mass balance stakes. The data coverage is 14 about 14% of the whole glacier area. The results from GPR data show a maximum depth of 235±18 m and a decreasing 15 mean depth ranging from 134±18 m in 2013 to 110±18 m in 2020. Additionally, we estimate a mean ice thickness change 16 rate of -4.2±3.2 m yr -1 between 2014 and 2020 with GPR data alone, which is in agreement with the elevation change in the 17 same period. The latter was estimated with the more accurate surface elevation data, yielding a change rate of -3.2±0.2 myr -1 , 18 and hence, confirming a negative glacier mass balance. The data set can be valuable for further analysis when combined with 19 other data types, and as input for glacier dynamics modeling, ice volume estimations, and GLOF (glacial lake outburst flood) 20 risk assessment. The complete dataset is available at https://doi.org/10.5281/zenodo.5571081 (Oberreuter et al, 2021). 21


Introduction 22
The glacier variations are sensitive to climate change (Oerlemans, 2005), particularly that of tropical glaciers (Rabatel et   The GPR data coverage from 2013 to 2020 is described in Figure 2 and in Table 1. Data from years 2013-2017 were 64 collected by the Autoridad Nacional de Agua (ANA), while data from years 2018-2020 were collected by Instituto Nacional 65 de Investigación en Glaciares y Ecosistemas de Montaña (INAIGEM). All data were obtained with the same radar system 66   two bistatic-shaped antennas. Its main features are detailed in Table 2. The transmission unit is able to generate a high-77 voltage impulse at a pulse repetition frequency (PRF) of 1 kHz, with an output of 1.4 kV. Its two integrated GPS/GLONASS 78 antennae (with GPS navigator accuracy) allow the trigger synchronization between receiver and transmitter, which is 79 important for setting the beginning of each trace. The two antennas are dipole Wu-King type (Wu & King, 1965), with 5 80 MHz of central frequency and bandwidth. During the capture process, the raw data are transferred to a handheld rugged 81 computer (PDA) for storing and real-time visualization of data. For a better deploying of the GPR, a group of five people 82 was needed during the surveys as shown in Figure 3. 83 The receiver unit digitalizes the signal at a sampling rate of 80 MHz with 16 bits of resolution, yielding a sample length (or 86 time increment) of 12.5 ns, and is able to stack from 256 to 4096 traces, with a trace length from 256 to 1024 samples. The 87 digital system enables the user to set a fixed trace length (in samples) between 256 and 1024 (2 n format). The proper choice 88 of the trace length depends on the estimated depth, ranging from 268 m to 1075 m, assuming a depth-averaged wave 89 propagation velocity in ice (c) of 0.168 m/ns (Glen & Paren, 1975). 90 With the interpretation of the two-wave travel time ( ), we estimate the ice depth ( ) using equation (1). 91 (1)

GPR data processing 95
The first stage of the processing considers data preview and format conversion with the software RADAR View, which is the 96 software provided by the GPR company. Then, the files are able to be imported into software Reflexw, where profiles are 97 loaded and integrated with GPS data, and processing takes place, including the following steps: 98 a) Normalization of amplitude-saturated traces: Due to different sources of interference, traces are sometimes affected by 99 saturation in the signal amplitude, which leads visually to discontinuities in the profile. In this case, the normalization 100 consists of applying the calculation Amplitude/max(Amplitude) to every trace. 101 b) Dewow filter: Low-pass filter used to eliminate the "wow" effect, caused by a low frequency component in the signal. 102 The filter is defined by a window length which should include the first arrival waveform. In this case, the windo w length 103 parameter was set to 100 ns, which satisfies the requirement. 104 c) Butterworth bandpass filter: To reduce the noise that is present in the signal outside the frequency range, where the radar 105 operates, which is centered at 5 MHz. The lower and higher boundaries of the filter were set to 2and 10 MHz, respectively. 106 d) Definition of zero time: The direct wave travelling through the air between source and receiver should be eliminated in 107 order to estimate correctly the two wave travel times of the corresponding reflections. In this case, the zero time was equal 108 for each dataset varying from 305 ns to 380 ns. 109 e) Trace interpolation: To have a better visualization of the profile and as a necessary step for the migration, an equidistant 110 trace interpolation was performed. Here a 1 meter separation between traces was used, preserving a similar number of traces 111 per profile, and consequently, keeping the same horizontal resolution. 112 f) Migration: In order to set the reflectors to their real position, a migration process was performed. It was used a 1D 113 Kirchoff migration because of its better signal to noise ratio in this particular case. Here the migration profile window was 114 set to 250 traces and the depth-average electromagnetic wave velocity was set to 0.168 m/ns. This velocity was in agreement 115 with hyperbolae shapes. 116 g) Topographic correction: This procedure enables to move the traces up and down in the profile according to the elevation 117 that is store in the trace header, which was defined when importing the GPS data to the profile The results enable to visualize 118 correctly the surface elevation and the bedrock elevation along the profile (see lower panel of Figure 4). 119 After processing, the bedrock interface was manually picked at each radar profile and then exported to GIS software. The 120 processed radargrams are also visualized in OpendTect in order to examine the crossover differences. An example of two 121 processed profiles is shown in Figure 4 Table 3. An automated 134 total station is a total station that automatically and quickly follows the target (prism) with a laser while the radial 135 topographic survey is being performed, facilitating the task in the fieldwork. The radial topographic survey is a source-to-136 receiver procedure, which aims to determine the position of several points over a surface, where the source is a total station 137 with known reference position (XYZ coordinate) and elevation over the ground and the receiver is a prism with known 138 https://doi.org/10.5194/essd-2021-336 elevation over the surface. The total station estimates the source-receiver distance and angles in order to calculate the 139 position of the receiver in relation to the source position. In this study, one single source reference position was used for each 140 year, which was located outside the glacier area, at a fixed rock position of coordinates E=209143.7, N= 9007853.1, datum 141 WGS84, UTM18S, elevation: 4764.1 m.a.s.l. The elevation of the total station was estimated with navigator GPS and it was 142 used the same of each year. Then, each source-receiver surveyed with prism is stored within the memory of the total station. 143 Once the radial topographic survey has been completed, the data is downloaded and processed . covering between 11% and 18% of the total glacier area (Table 3).When measuring the distances between nearest points of 149 every dataset, a mean value of ~25 m was found. In order to obtain the DEMs, the topographic points were converted to a 150 5m by 5 m raster and interpolated using the Inverse Distance Weighted interpolation tool in ArcGIS. 151 Due to better resolution and timespan of data, the first dataset was used to estimate the surface elevation changes of the 152 glacier, by subtracting the digital elevation model of 2020 from that of 2012. 153 In order to take advantage of the radar data coverage, the bedrock elevation was estimated by subtracting the ice thickness 154 interpretation from the surface elevation (base on the GPS of the radar survey) in every profile of the GPR data. 155

Errors 157
The estimation of errors in ice thickness (ε Hdata ) was obtained following Lapazaran et al. (2016), who split it into two 158 components: a) the error in the value of ice thickness due to GPR measurement (ε HGP R ), without taking into account where it 159 was obtained and b) the error in ice thickness due to uncertainties in horizontal positioning (ε Hxy ). ε Hdata is variable and 160 hence, different for every GPR measurement point. 161 The term ε HGP R is a function of the propagation velocity in ice (c) and its uncertainty (ε c ), TWTT, and timing error (ε τ ).

162
The velocity c is dependent of water content and ice purity. In this case, we set it to 0.168 m/ns as used by Gacitúa

165
The term ε τ is a function of the central frequency of the GPR (see Table 2) and c, which adds up 16.8 m.

166
The term ε Hxy should be calculated considering the precision of GPS, speed of transportation and geometric aspects 167 (Lapazaran et al, 2016). ε Hxy is then calculated a posteriori, unlike the other errors mentioned in this section, because an ice 168 thickness measurement is needed in order to estimate the error that applies on that measurement. In this case, we considered 169 5 meters of uncertainty in horizontal positioning because of the use of single frequency GPS (there was no dual frequency 170 GPS). With this uncertainty in the XY plane, a 5m by 5m cell point -to-raster for ice thickness procedure was performed, 171 estimating the difference between maximum and minimum interpreted ice thickness within each grid cell and then defining 172 this value as ε Hxy /√ .

173
An example of the estimated errors in ice thickness ε HGP R , ε Hxy and ε Hdata for 2014 dataset is shown in Figure 5. Also, a 174 crossover analysis for the same dataset and for the whole ice thickness dataset is provided in Table 4. The analysis shows a 175 maximu m difference of 10 meters, which is less than the radar resolution ε HGP R.

Surface elevation 193
The glacier surface elevations along profile CC' (Figure 1), which were obtained from interpolated DEMs and mass balance 194 stakes between 2012 and 2020, are presented in Figure 6(a). The elevation from available stakes measurements are in 195 agreement with that from DEMs (both data types are found in years 2014 and 2018). 196 In terms of elevation change for the periods 2012-2014, 2014-2018 and 2018-2020, mean differences of -5.8, -12.0 and -7.9 199 m were found, in the common area of the DEMs which mainly includes the ablation area. This leads to thinning rates of -200 2.9±1.0 m yr -1 , -3.0±0.9 m yr -1 , and -4.0±1.8 m yr -1 , respectively (Table 5). 201 On the other hand, for the whole study period 2012-2020, a mean difference of -25.6 m with a standard deviation of 3.2 m 202 (13% of the absolute mean value) was found, as shown in Figure 6(b) and in Table 5. Assuming a temporal error of 0.

Ice thickness data 210
The interpreted ice thickness data at Artesonraju glacier from 2013 to 2020 as well as the 2014-2020 thickness difference are 211 presented in Figure 7, while the main statistics of the interpreted data can be found in Table 6.

Conclusions 246
We provide a robust dataset contributing with interpreted ice thickness data and ice surface elevation data at the ablation area 247 of Artesonraju glacier in Caraz-Áncash, from 2012 to 2020. Data analysis reveals two aspects: a) a strong decadal mass 248 balance reduction and b) the existence of an overdeepening up to 700 m upstream the central flowline, which may lead to an 249 expansion of lake Artensoncocha Alta up to that point. Data here provided is useful as input for glacier dynamics modelling 250 and can be complemented with other data types for further analysis. It is also necess ary to assess the risks associated with the 251 volume increase of the proglacial lake Artesoncocha Alta and with GLOF generation, which may constitute a hazard for the 252 communities located downstream the Artesonraju glacier.