Probabilistic estimation of glacier volume and glacier bed topography: the Andean glacier Huayna West

Glacier retreat will increase sea level and decrease fresh water availability. Glacier retreat will also induce morphologic and hydrologic changes due to the formation of glacial lakes. Hence, it is important not only to estimate glacier volume, but also to understand the spatial distribution of ice thickness. There are several approaches for 5 estimating glacier volume and glacier thickness. However, it is not possible to select an optimal approach that works for all locations. It is important to analyse the relation between the different glacier volume estimations and to provide confidence intervals of a given solution. The present study presents a probabilistic approach for estimating glacier volume and its confidence interval. Glacier volume of the Andean glacier 10 Huayna West was estimated according to different scaling relations. Besides, glacier volume and glacier thickness were estimated assuming plastic behaviour. The present study also analysed the influence of considering a variable glacier density due to ice firn densification. It was found that the different estimations are described by a lognormal probability distribution. Considering a confidence level of 90 %, the estimated 15 glacier volume is 0.0275km ± 0.0052km. Considering a confidence level of 90 %, the estimated glacier thickness is 24.98 m with a confidence of ±4.67 m. The mean basal shear stress considering plastic behaviour is 82.5 kPa. The reconstruction of glacier bed topography showed the future formation of a glacier lake with a maximum depth of 32 m. 20

estimating glacier volume and glacier thickness. However, it is not possible to select an optimal approach that works for all locations. It is important to analyse the relation between the different glacier volume estimations and to provide confidence intervals of a given solution. The present study presents a probabilistic approach for estimating glacier volume and its confidence interval. Glacier volume of the Andean glacier 10 Huayna West was estimated according to different scaling relations. Besides, glacier volume and glacier thickness were estimated assuming plastic behaviour. The present study also analysed the influence of considering a variable glacier density due to ice firn densification. It was found that the different estimations are described by a lognormal probability distribution. Considering a confidence level of 90 %, the estimated

Introduction
Glaciers may be considered as the most important water reservoirs since they store 68.7 % of the world's total fresh water (Gleick, 1996). Unfortunately, they are retreating. Glacier area loss was observed in the Alps , in the Andes (Ramirez et al., 2001), in Africa (Kaser et al., 2004), in the Himalayas (Naito et al., 2012), in 25 Oceania (Hoelze et al., 2007), in the Antarctic (Rignot, 1998) and in the Arctic (Lenaerts et al., 2013). Different studies showed that glacier retreat will influence sea level rise and water shortage. Dyurgeyov and Meier (2005) estimated that mountain glaciers may induce a sea level rise of approximately 0.65 m. This level rise may result in a spatial shift of coastal lines (Crooks, 2004) and flooding of some delta areas (Dissanayake et al., 2009). Glacier retreat will also affect water resources. Glacier retreat will induce 5 changes in the streamflow that will affect the water availability downstream (Vuille et al., 2008). Thus, it is important to estimate glacier volume in order to predict future water availability and sea level rise (Kaser et al., 2010;Raper and Braithwaite, 2006).
Field measurement campaigns are a direct method for estimating glacier volume. Glacier thickness and glacier volume can be estimated using ground-penetrating radar 10 (GPR) (Singh et al., 2010;Degenhardt, 2009;Monnier et al., 2011;Navarro et al., 2005;Lee et al., 2010) and radio echo sounding (RES) (Andreassen et al., 2012;Jacobell et al., 2010;Dowdeswell and Evans, 2004;Goodman, 1975). However, GPR and RES are impractical methods for remote and extensive glaciers. Hence, new and simpler alternatives for estimating glacier thickness and glacier volume have been de- 15 veloped. Farinotti et al. (2009) combined mass balance and ice-flow dynamics. Mathematical ice-flow models were developed as an analytical alternative that allows estimating glacier thickness (Colgan et al., 2012;Goldberg and Sergienko, 2011). Such methods require glacier velocity, which is also difficult to obtain; in fact, some studies use ice-flow models for the estimation of glacier velocity (Wang et al., 2012;Debol-20 skaya and Isaenov, 2010; Ren et al., 2012;Calvo et al., 2010). Other studies combine such models with ground-measured data (Farinotti et al., 2013;Zekollari et al., 2013;Colgan et al., 2012). Other studies use analytical models for the estimation of mass balance (Michel et al., 2013;Morlighem et al., 2011). Kääb (2000) combined an iceflow model with photogrammetric data for the reconstruction of glacier mass balance 25 in the Swiss Alps. A different approach is the combination of field-measured data with regression techniques in order to develop empirical relations between glacier thickness and certain local characteristics (Peduzzi et al., 2010;Clarke et al., 2009). However, the latest method lacks a physical background and can only be applied locally.

3933
Other studies tried to develop universal methods. Bahr et al. (1997) developed a physically based volume-area power-law scaling relation. This relation is based on ice dynamics constraints due to ice rheology and typical climatic-topographic conditions of glacierized areas. Nevertheless, such a relation depends on two empirical parameters. Several studies suggested different values for such parameters (Chen and 5 Ohmura, 1990;Meier and Bahr, 1996;Baraer et al., 2012). Thus, this method is prone to an important degree of uncertainty. Nevertheless, this is the most practical method for estimating glacier volume. Besides, some studies suggest that the initial uncertainties of scaling methods are small when analysing future glacier evolution (Radic et al., 2008). One limitation of this method is that it only predicts the total glacier volume, neglecting the distributed thickness.
Spatially distributed glacier thickness is valuable information since it would allow us to predict the glacier bed topography (GBT) (Binder et al., 2009). Such a GBT is important for modelling the glacier evolution (Huss et al., 2010), for estimating possible changes in the runoff regime (Huss et al., 2008), and for predicting the formation of 15 future lakes (Frey et al., 2010). According to Uchupi et al. (2001), the morphology imparted by glaciation has important consequences on drainage characteristics and in some cases may regulate the discharge. Some studies suggest that sub-glacial processes are related to volcanic activity (Scharrer et al., 2008). Other studies suggest that sub-glacial drainage may have storage and delayed discharge, influencing the 20 flow of the ice sheet (Baelum and Benn, 2011;Thoma et al., 2010). Some glacial lakes may be covered by glaciers, but after glacier retreat such lakes are exposed . Glacial lakes may overflow and fail, causing outburst floods (Tadono et al., 2012). Actually, several flood events due to glacial lakes have already been reported .
A new and more practical approach for estimating glacier ice thickness and glacier volume is the glacier bed topography (GlabTop) approach proposed by Paul and Linsbauer (2012). Assuming a plastic behaviour, glaciers flow easily enough to redistribute mass and prevent stresses from rising above a given limit (Cuffey and Paterson, 2010; Nye, 1967); thus, the basal shear stress is equivalent to the plastic yield stress. Such a plastic assumption is supported by field measurements showing that glacier deformation is best reproduced considering a plastic deformation (Kavanaugh and Clarke, 2006). GlabTop assumes the glacier thickness as a function of surface slope and basal shear stress. Such assumptions were successfully applied to glaciers in different parts 5 of the world (Li et al., 2012;Paul and Linsbauer, 2012;Linsbauer et al., 2012). Recently, Clarke et al. (2013) proposed to estimate a basal shear stress for different glacier flow units defined by their ice catchment, and then to estimate glacier thickness. However, the latest method requires the ice flux, which is a difficult value to obtain.
In the present study we assume that there is not a unique absolute solution for es-10 timating glacier volume; all the proposed methods have some degree of probability and neglecting some estimations may lead to biased conclusions. Glacier volume and glacier thickness of the Andean glacier Huayna West were estimated by the different proposed scaling relations. Besides, distributed glacier thickness and glacier volume were estimated by the GlabTop method, considering different possible values of 15 basal shear stress and analysing the differences between the consideration of a constant glacier density and depth-variable glacier density. Then, the estimations were statistically compared and analysed. It was found that the different volume estimation methods are best explained by a lognormal probability distribution function. Then, the probable glacier volume and glacier thickness were expressed considering appropriate 20 confidence intervals. With the spatially distributed glacier thickness, it was possible to reconstruct the glacier bed topography (GBT). The GBT shows the formation of a future lake. The estimated area-depth ratio of such a lake fits reasonably with observed ratios of glacial lakes from other locations.

Study area 25
The study area is the Huayna West glacier in the Bolivian Andes. It is located on the west side of Huayna Potosí, a massif mountain with a maximum elevation of 6088 m 3935 above sea level (m a.s.l.). Huayna Potosí is located 30 km north of La Paz, Bolivia. It separates the dry Altiplano in the west from the wet Amazonian Basin in the east. Due to these climatological differences, the equilibrium line altitude (ELA) of the west side (5700 m a.s.l.) is shifted almost 1 km above the eastern ELA (Condom et al., 2007). Huayna glacier is located in the Capricorn tropic, where the climate is characterized by 5 two seasons with a period of precipitation and a dry period (Mote and Kaser, 2007). It has a wet season during the austral summer (November-March) and a dry season during the austral winter (April-October). This glacier is highly related to the water supply of La Paz-El Alto conurbation (Bolivia), and is currently being studied within the Glacier Retreat impact Assessment and National policy DEvelopment (GRANDE) 10 project.

Methodology
Glaciers were delineated using remote sensing data from the Advanced Visible and Near Infrared Radiometer type 2 sensor (AVNIR-2) of the Advanced Land Observation Satellite (ALOS). The observation date of the image was 6 May 2010. The ALOS 15 AVNIR-2 images were analysed and processed with the multispectral image data analysis system software (Landgrebe, 2005;Dundar and Landgrebe, 2004). The images were classified using an unsupervised classification that allowed distinguishing into four different coverage types: impoundments, rocky ground, bare soil and glaciers ( Fig. 1).
It is important to note that the ALOS AVNIR-2 image is more recent than the respective 20 image used for the Randolph Glacier Inventory 3.0, which was observed on 31 May 2003. Besides, it was not possible to use the GLIMS database, since such a database does not include this glacier (Raup et al., 2007). The flow lines were obtained by processing the global digital elevation model (GDEM) provided by the advanced spaceborne thermal emission and reflection ra-25 diometer (ASTER). It is important to remember that surface slope is the dominant factor on bed stress (Nye, 1954). Thus, the surface slope α may be assumed equal to the bed slope (Clarke et al., 2013). Although some studies used SRTM (Shuttle Radar Topography Mission) data for studying glaciers and estimating glacier volume (Cooper et al., 2010;Surazakov and Aizen, 2006), in the present study we used the GDEM ASTER since it is a more recent product with better resolution and provides an accurate delineation of the study area. The DEM was processed with the TauDEM 5 algorithm implemented in the GIS software MAPWINDOWS (Tarboton, 1997). Overlaying the delineated basin over a remote sensing image of the study area shows that the delineated basin accurately matches the topography (Fig. 2).
Then, glacier thickness was estimated at the flow-path lines assuming perfect plasticity and the GlabTop approach (Linsbauer et al., 2012) described by Eq. (1).
where h is the glacier thickness (m), τ the basal shear stress (kPa), g the gravity acceleration (9.79 m s −2 ), and θ the slope ( • ). The most popular estimation of basal shear stress is to consider it as a function of the elevation range of the glacier (Eqs. 2 and 3); however, it is important to note that since tropical glaciers have high mass 15 balance gradients, it is reasonable to expect a higher basal shear stress described by Eq. (4), instead of Eq.
(3) (Haeberli and Hoelze, 1995). 20 where ∆H (km) is the elevation range of the glacier. Thus, the basal shear stress may be any value between the limits imposed by the above-mentioned equations. In the present study we estimated glacier depth considering both limits and the average. The first basal shear stress from Eq. (2) was denoted 25 as τ1 (64.20 kPa), the second basal shear stress from Eq. (4) as τ2 (136.50 kPa), and the average basal shear stress was denoted as τ3 (100.35 kPa). 3937 Glacier density is usually assumed as a constant value of 900 (kg m −3 ) independent of the depth (Linsbauer et al., 2012;Paul and Linsbauer, 2012;Li et al., 2012). However, it is important to consider that the progressive transformation of snow into ice describes a density-depth relation (Cuffey and Patterson, 2010). Recent studies state that the assumption of a constant density may induce important errors in geodetic 5 studies (Huss, 2013). Hence, in the present study both cases were considered: a constant density and a depth-variable density. The density-depth relation was assumed to respond to a parabolic equation (Shumsky, 1960). The density-depth relation was obtained using density measurements from different Andean glaciers (Ginot, 2001). The density-depth relation was described by Eq. (5).
where d is the glacier depth (m) and ρ is the glacier density (kg m −3 ) at such depth. However, such a relation provides the density at a specific depth. In order to get an average density of a glacier column of a given depth, an additional relation was obtained considering that any computed density at a given depth remains constant for intervals 15 of 1 m. Average density for different depths was then computed by where ρ k is the average density for a glacier column of depth k , ρ i the glacier density at depth i and i is an integer number. Then, the glacier thickness points were geo-referenced and the distributed glacier 20 thickness was estimated by applying a kriging interpolation, which is an interpolation technique widely used in glaciological and hydrological studies (Bamber et al., 2013(Bamber et al., , 2009Binaghi et al., 2013;Sørensen et al., 2011). Before the interpolation, the outline of the glacier was assumed to have a thickness of 1 m. This assumption was done in order to limit the boundaries to be interpolated and to avoid possible negative values.
The distributed glacier thickness maps allowed obtaining the glacier volume and average thickness. The glacier volume was obtained by multiplying each thickness by its area (900 m 2 ). The average thickness was estimated by averaging the glacier thickness values. Considering the two possible densities and the three possible values of the basal shear stress, we have six possible glacier volumes and six possible average 5 thickness values. In order to evaluate the results, glacier volume was estimated by other available methods. The most popular method is the area-volume scaling relationship . This method relates glacier volume and glacier area by a power-law scaling relation (Eq. 7).
where A is the area of the glacier (km 2 ), V the volume of the glacier (km 3 ), γ the scaling exponent and c a proportional constant. The value of gamma depends on the geometry of the glaciers. Different studied glaciers from different locations developed volume-area scaling relations with different values between 0.02 and 0.0597 for c and 15 between 1.12 and 1.5 for γ. (Macheret et al., 1984;Yafeng et al., 1981;Meier and Bahr, 1996;Radic and Hock, 2010;Grinsted, 2013). Adhikari and Marshall (2012) suggest different equations according to the glacier mean slope, glacier area and shape factor defined as the width-length ratio. Besides, recent studies showed that γ also depends on the transient state of the glaciers; glaciers under warmer scenarios may have γ 20 values higher than 2.0 (Radic et al., 2007). Analysing area volume data of the last 60 yr of the extinct glacier Chacaltaya (Francou et al., 2000;Ramirez et al., 2001), we found a γ value of 2.0207 and a c value of 0.1091. Recent studies stressed the importance of calculating glacier volume by considering thickness-area scaling relationships (Eq. 8). Thus, in the present study we also con-25 sidered four thickness-area relationships (Huss and Farinotti, 2012;Ohara et al., 2013;Bodin et al., 2010;Nicholson et al., 2009). Three of such relationships were previously 3939 developed and applied to South American Andean glaciers.
whereh is the mean glacier thickness (m). Then, multiplying the mean thickness by the area gives the volume. Considering the 6 initial possible volumes, the 24 volumes from the volume-area 5 relationships and the 4 volumes from the thickness-area relationships, we have a total of 34 possible glacier volumes (Table 1).
Since the selection of an arbitrary single solution may lead to bias and wrong estimations, a statistical analysis was performed in order to get the most probable solution. Thus, a statistical analysis was performed in order to find the probability distribution 10 function (PDF) that describes the best of the different volume estimations. Seven probability distribution functions were considered: -Beta distribution is described by Eq. (9).
where α 1 and α 2 are shape parameters of the distribution, a and b boundary 15 parameters, B the beta function and x the variable.
where γ is the scale parameter of the distribution.
where α and β are the shape parameters of the distribution, and Γ is the gamma function.
where σ is the scale parameter and µ the location parameter of the distribution.
The distributions were evaluated considering the tests of chi-square, Kolmogorov-Smirnov and Anderson-Darling: the chi-square test compares the observed frequencies with the frequencies of an assumed theoretical distribution (Ang and Tang, 1975).
The Kolmogorov-Smirnov test performs a comparison between the experimental cu-15 mulative frequency and an assumed theoretical distribution (Marsaglia et al., 2003). 3941 The Anderson-Darling test also compares the cumulative frequency, but it gives more weight to the tails (Anderson and Darling, 1954). The statistical fit analysis was performed with the software EasyFit 5.5 (http://www.mathwave.com/easyfit-distributionfitting.html). Once the best fit distribution was selected, it was possible to associate probabilities t statistic (Brownlee, 1960). However, such a method provides good confidence intervals for normally distributed population; the method may lead to wrong predictions for a population with a different distribution (Singh and Nocerino, 2002). Several studies analysed different methods for estimating confidence intervals for data with different probability distribution functions (Endo et al., 2009;Gutierrez et al., 2007;Wu et al., 15 2005; Takada and Nagata, 1995). The method used for the estimation of confidence interval was selected according to the results and will be explained in the respective section.
The glacier volume was estimated as the mean value of such a distribution. Then, a trial and error approach was used in order to find the basal shear stress that provides 20 a volume that fits the estimated mean glacier volume. Then, the estimated glacier thickness (GT) was subtracted from the glacier surface elevation (GSE) in order to get the glacier bed topography elevation (GBTE): where the subscripts i and j identify the glacier cells according to their row and column 25 in the raster. GSE was obtained from the DEM.

Results and discussion
The range of possible basal shear stress values is between 64.20 kPa and 136.50 kPa. The value of 64.20 kPa was denoted as basal shear stress 1 (τ1), the value of 136.50 kPa as basal shear stress 2 (τ2), and the average value of 100.35 kPa as basal shear stress 3 (τ3). Figure 3 shows a comparison of the distributed glacier thickness 5 considering constant density and variable density, and the different basal shear stress considered. The value of the basal shear stress has a strong influence on the total volume estimation. The results using a basal shear stress of 64.20 kPa are about half the value of the results using a basal shear stress of 136.50 kPa. Such results are not unexpected ones since τ1 is almost 50 % of τ2, and the basal shear stress is directly pro-10 portional to the glacier thickness. The differences between the two density approaches are much lower. The variable density provides higher volume estimations. Such differences are expected since the variable density is lower than 900 kg m −3 . Since glacier density is in the divisor of Eq. (1), lower densities will provide higher results. Considering variable density and τ1, the total volume is 10.5 % higher than considering constant 15 density and τ1. Considering variable density and τ2, the total volume is 3.3 % higher than considering constant density and τ2. Considering variable density and τ3, the total volume is 3.9 % higher than considering constant density and τ3. However, the uncertainties from glacier density are minor compared with the uncertainties from the basal shear stress. The biggest differences between a constant and a variable glacier density 20 are in the spatial distribution of the glacier thickness. Figure 3 also shows that the area covered by low thickness is bigger when considering constant density. In all the cases the deepest part of the glacier is located at some 180 m from the east boundary. The deepest part is elongated with a north northeast-south southwest direction and a total longitude of 370 m.

25
The 34 methods provide different volume estimations ( Table 2). The lowest estimation is provided by Bahr et al. (1997), and the highest is provided using the coefficients deduced from Francou et al. (2000). The average estimation is 0.027, which is close to 3943 the estimations according to Shiyin et al. (2003) and Macheret et al. (1984). However, such a mean value considers a simple normal average. After fitting the volume estimation to the different probability distribution functions (Table 3), it was found that the lognormal distribution is the one that has the best performance under the three best-fit tests ( Table 4). The gamma distribution is the one with the second-best performance. 5 The distributions that ranked as the third and the fourth are the logistic distribution and the Weibull distribution. The normal distribution ranks as the fifth distribution. The worst performance is from the beta distribution and the exponential distribution. Table 5 shows the mean and standard deviation (std) values of the different PDFs. A superficial analysis of the result may lead to the idea that the differences are quite low.
If we consider the lognormal mean value as the target value, the mean estimation from the beta distribution has the highest difference with an overestimation on 12.60 %. The Weibull distribution provides an underestimation of 3.96 %, and the other distributions have an overestimation about 0.25 %. However, the standard deviation shows much higher differences. The beta distribution std is 37.28 % higher, the Weibull std 15.78 % 15 lower, the exponential std almost three times higher, and the std from the other distributions 9.45 % higher. Such differences have a strong influence over the confidence intervals. Considering the lognormal distribution the Huayna glacier has an estimated volume of 0.0275 km 3 . Such a volume is equivalent to a mean thickness of 24.98 m. This volume can be obtained by considering a basal shear stress of 82.5 kPa. 20 It is important to point out that the US Environmental Protection Agency (EPA) performed an evaluation of different methods for the estimation of confidence intervals for data with different distributions (i.e. not normal distributed data) and suggested the use of different methods such as the Chebyshev approach (Singh et al., 2004). Actually, other studies also suggest the use of the Chebyshev approach for the estimation of confidence intervals of other distributions and unknown distributions (Amidan et al., 2005). The Chebyshev theorem tells us that, for any data set, the proportion of data (pd) that lies within k standard deviation on either side of the mean is at least (Almukka-hal et al., 2011) Then, the upper confidence interval (UCI) and lower confidence interval (LCI) were obtained by adding and subtracting the respective confidence.
where n is the number of data values. Table 6 shows the confidence values (Conf) to be added and subtracted in order to get the confidence interval for different probabilities.
Using the proper probability density function allows to obtain the cumulative distribution function (CDF) (Fig. 4), useful for estimating the probability that the variable, in this case the volume, is at least a given value. Models 11, 10, 13 and 6 are within the 10 % CDF; thus, there is a 90 % probability that the glacier volume is higher than those values. Models 27 and 1 are in the 20 % CDF. Models 12,33,9,28,8  outlier. Figure 5 shows the reconstructed glacier bed topography (GBT). The GBT shows the formation of a future lake in the quadrant D-3 once the glacier disappears. This lake has an area of 0.07 km 2 and a maximum depth estimated of 32 m. The areamaximum depth relation of the glacier is a reasonable value that fits reasonably with 25 3945 other estimations. For instance, Sakai (2012) developed a power area-depth relation considering several glacial lakes. Applying such a relation to the present lake gives a maximum depth of 26.06 m.

Conclusions
Theoretical approaches for glacier volume estimation are influenced by coefficients that 5 depend on local conditions. They cannot predict with certainty the volume of a given glacier; however, they can predict the glacier volume with an associated probability. This study presented a comparison of different volume estimations and a statistical approach for predicting the associated probability of the estimated glacier volume. Besides, the influence of glacier density was evaluated. 10 The basal shear stress has a strong influence in estimation of the glacier thickness and glacier volume. The influence of the glacier depth is small compared with the influence of the basal shear stress influence for estimating the total glacier volume. Nevertheless, the glacier density has an important influence in the spatial distribution of the glacier thickness. Assuming a constant glacier density, the predicted glacier thick-15 ness is lower than the predicted assuming a variable glacier density, especially in lower areas.
The different glacier volume estimations are related by a lognormal probability distribution. The mean glacier volume is 0.0275 km 3 . Considering a 90 % confidence level, such an estimation has an uncertainty of 0.0052 km 3 . The Huayna West glacier has 20 a mean thickness of 24.98 m. Considering a 90 % confidence level, such an estimation has an uncertainty of 4.67 m. The mean glacier volume could be obtained assuming plastic behaviour and a basal shear stress of 82.50 kPa. This is a reasonable result since it is within the possible ranges of basal shear stress. Such basal shear stress may be considered as represen-25 tative of the glacier.
The equations proposed by Macheret et al. (1984) and Shiyin et al. (2003) provide the nearest estimation to the mean glacier volume. The glacier volume estimated with the GlabTop approach considering basal shear stress of 64.20 kPa and variable glacier density is the closest to the lower confidence interval. This estimation has a cumulative distribution function (CDF) of 35 %. The glacier volume estimated with the GlabTop ap-5 proach considering a basal shear stress of 100.35 kPa and variable glacier density is the closest to the upper confidence interval. This estimation has a cumulative distribution function (CDF) of 73 %.
Glacier bed topography showed the formation of a future glacial lake. The estimated area and depth of this lake have a reasonable agreement with dimensions observed at 10 other glacial lakes.
Glacier retreat will not only influence the water availability, but also induce morphometric changes to the hydrological drainage network.