The Cryosphere Assessing high altitude glacier thickness , volume and area changes using field , GIS and remote sensing techniques : the case of Nevado Coropuna ( Peru )

Higher temperatures and changes in precipitation patterns have induced an acute decrease in Andean glaciers, thus leading to additional stress on water supply. To adapt to climate changes, local governments need information on the rate of glacier area and volume losses and on current ice thickness. Remote sensing analyses of Coropuna glacier (Peru) delineate an acute glaciated area decline between 1955 and 2008. We tested how volume changes can be estimated with remote sensing and GIS techniques using digital elevation models derived from both topographic maps and satellite images. Ice thickness was measured in 2004 using a Ground Penetrating Radar coupled with a Ground Positioning System during a field expedition. It provided profiles of ice thickness on different slopes, orientations and altitudes. These were used to model the current glacier volume using Geographical Information System and statistical multiple regression techniques. The results revealed a significant glacier volume loss; however the uncertainty is higher than the measured volume loss. We also provided an estimate of the remaining volume. The field study provided the scientific evidence needed by COPASA, a local Peruvian NGO, and GTZ, the German international cooperation agency, in order to alert local governments and communities and guide them in adopting new climate change adaptation policies. Correspondence to: P. Peduzzi (pascal.peduzzi@unepgrid.ch)


General context
Changes in glaciers and ice caps are good indicators of climate change (Zemp et al., 2008) and the current trend shows that a majority of the world glaciers have undergone a reduction in their mass at an accelerating rate.The mass loss was greater in the period 1990/91 to 2003/04 than in the previous period 1960/61 to 1989/90 (Bates et al., 2008).This is of concern given that about one-sixth of the world's population depend on glacier and snow melting for their water supply (Bradley et al., 2006).
In Peru, the population growth and rising water demand for agriculture, domestic and economic activities generate an increased pressure on water resources.As the rainy season is concentrated during four months of the year, the role of glaciers is crucial for spreading out the water supply during the dry season.Higher limit between rain and snow precipitation reduces the buffering role of ice and snow, thus increasing flood risk during the wet season and reducing dryseason water supplies.This is of concern particularly in China, India and Asia, but also in the South American Andes, where a large fraction of the population relies on the glacial melt for water supply and hydropower (Barnett et al., 2005).In the South American region, the glacier monitoring for the period 1970-1996 revealed an acute retreat of Andean glaciers, with glacier coverage decreasing from 725 km 2 in 1970 to 60 km 2 in 1996 in Cordillera Blanca, Peru (Silverio and Jaquet, 2005).

Assessing change of ice volume in Nevado
Coropuna, Peru (6500 m a.s.l.) The present study on the Coropuna Glacier was made at the request of the Cooperación Peruana Alemana de Seguridad Alimentaria (COPASA) Special Project of Arequipa Regional Government, in collaboration with the Deutsche Gesellschaft für Technische Zusammenarbeit (GTZ).They asked for scientific-based evidence of glacier area and volume changes, in order to assess whether climate change adaptation policies on water supply should be introduced at local government and communities levels.
The study was carried out by a team from UNEP/GRID-Europe and the University of Geneva.It assessed glacial retreat using both satellite imagery analysis and in situ measurements of the Coropuna Glacier.
This paper describes how glacier area was monitored through time and how to measure the change in volume of the glacier using different Digital Elevation Models (DEMs) as well as evaluate the current (2004) ice thickness.To this end, a field mission was carried out using a Ground Penetrating Radar (GPR) coupled with a Ground Positioning System (GPS).It provided profiles of ice thickness on different slopes, orientations and altitudes.These profiles were used for modelling the entire remaining glacier volume using Geographical Information System (GIS) and statistical multiple regression techniques.
Given the limited financial resources of the local governments and development organisations in Peru, simple and low-cost techniques to measure changes in glacier volume and the remaining ice volume were needed.The purpose of this paper is to test how, with limited budgets and using locally available hardware, scientific evidence of glacier area and volume variation as well as ice thickness can be obtained.2 Study area and data sources

Study area
The Coropuna Glacier is the third highest summit in Peru, culminating at 6446 m.It is located at 15.546 • S, 72.660 • W, about 155 km northwest of the city of Arequipa (Fig. 1).According to COPASA staff, 8000 people depend on the Coropuna Glacier for their water supply and it is estimated that 30 000 people depend indirectly on the glacier for their livelihoods.

Passive satellite sensors images
To compute the difference in glaciers area loss, we used a topographic map for the area in 1955 and four images from 1980, 1996, 2003 and 2008 (see Table 1).Only cloud-free images were used.Except the image 2003, which is still early in the season (May) and might still have some snow, all the others are taken far from the snow season.1955, 1997, 2000 and 2002.The DEM 1955 was generated by manually digitising the elevation contour lines from the topographic map of 1955.The DEM 1997 was purchased from the company SARMAP which produced it based on a pair of ERS-1 Synthetic Aperture Radar (SAR), a satellite from the European Space Agency.The DEM 2000 was based on radar measurements from the NASA Shuttle Radar Topographic Mission (SRTM, version 3).The DEM 2002 was derived from ASTER satellite data, DEM provided by United States Geological Survey (USGS) (see Table 2).

Field measurements
The purpose of the expedition was to measure the depth of the ice as well as taking GPS points for the adjustment of the DEMs.The 14 day-mission was undertaken between 13 and 26 August 2004.The team was composed of two scientists and 11 support staff.
The scientific instruments were chosen based on local availability (as opposed to most efficient).The GPR Ramac X3M included a 100 MHz shielded antennas.Much lighter GPR exist; however, it was the only GPR available in Arequipa.Three Global Positioning System receptors (GPS) and two regular office laptops for controlling the GPR unit and recording the data were used.Due to the limitations of the computer's hard disk at low pressure conditions (the reading heads would touch the disks and damage them), hard disks were removed.Computers and software were booted on CDs, and measurements were recorded on USB cards.

Estimation of ice thickness
Measuring the ice thickness was achieved using the GPR, with a sampling frequency of 438 MHz.Technical settings are specified in Table 3.The signal was assumed to travel through ice at 0.16 m/ns ± 0.01 m/ns according to other studies performed in similar conditions (Gruber, Ludwig and Moore, 1996;Moorman and Michel, 2000;Descloitres et al., 1999).The depth of ice can be measured by recording the time lag between the emission and the reception of the signal (see Eq. 1) Each recorded depth was coupled with geographical coordinates obtained by a GPS so that the profiles could be georeferenced and speed of radar over ground computed (given that for each GPS location, the time is also recorded).Due to time and access constraints, it was not possible to achieve a comprehensive coverage of the glaciated area.Instead, the mission proceeded along transects (shown on Fig. 2).The selected transects were chosen to provide samples including different altitudes, slopes, and aspects (slope orientation).We proposed (and tested) the hypothesis that these three variables would explain most of the ice thickness.Using multiple regression analysis, depth was modelled in areas where no measurements were taken.
During the mission, 10.6 km of transects were obtained from the GPR. Figure 2 shows the radar profiles.A GPR signal was recorded every two seconds, and a GPS location with hours, minutes and seconds approximately every 2 min.It was then possible to link GPS point with profile traces and accurately geo-reference the profiles.GPR transects were processed using the software "Reflex", "Ground Vision" and  .
"Kingdom Suite" to estimate ice depth.Due to the computer configuration that limited recording time windows, the bedrock was sometimes too deep to be detected (typically in volcanic craters, see Fig. 9).However, as an approximation, profiles can be extrapolated by following ice stratification.

Estimation of ice volume loss
Although passive satellite sensors, such as Landsat TM, provide an estimate of the area covered by ice (see Sect. 4.1) passive sensors do not provide information on the changes in ice thickness.However, the difference of ice volume may be estimated by using DEM time series.

Adjustment and corrections
In order to compare the different DEMs, several operations were needed.Firstly, all the DEMs were re-sampled to 30 m to compensate for different spatial resolutions.They were then reprojected in Universal Transverse Mercator (UTM), projection 18 south, datum WGS84.Finally, they were georeferenced so that they could be overlaid.This was performed using 32 control points, such as summits located outside the glacier area (on bare rocks).
Previous studies highlighted the issue of DEMs such as ASTER or SRTM which appear to have a bias following altitude (Berthier et al., 2004) or reveal significant errors when applied to rugged terrain and steep slopes (Kääb et al., 2002).In all DEMs, a reference area was chosen on surfaces that were not covered by ice, snow or vegetation (i.e.rocky or bare ground).This corresponds to the area outside the glaciated area of 1955.A sample of 1009 points located on an ellipse outside the glaciated areas was used to extract information on easting (X), northing (Y ), elevation, slope and aspect (slope orientation).These factors proved to influence the DEM errors in other studies (Gorokhovich and Voustianiouk, 2006).To verify whether the DEMs used were welladjusted, the differences in altitude were computed against the reference area (DEM 2000-DEM 1955;;DEM 1997-DEM 1955and DEM 2002-DEM 1955).
Then the differences between DEMs were plotted in 3-D along with X and Y .Despite the fact that the DEMs are all in the same geographic projection, Fig. 3  This can be corrected using the following linear regression (see Eq. 2): Where DEMt is the new values corrected for the X and Y .The weights a,b,c used to correct the DEMs are provided in Table 4.After this first correction, there were still some bias induced by aspect and elevation.To correct these, a quadratic equation was applied.The aspect variable was transformed by taking the cosine of the angle (expressed in radians).This is necessary as an orientation of 359 • is close to an orientation of 1 • .This distortion from elevation and aspect was corrected using Eq. ( 3): Where: DEMt2 = DEMt further corrected with elevations and cosine of the aspects El = Elevation As = Aspect (slope orientation) This second correction really improved the DEM 2002 (see Table 5), while for the DEM 1997 and the DEM 2000, the standard deviation increased, so the second correction was not applied to them.The weights used to obtain the DEM2002t2 were as followed: a = 0.1455, b = 0.0012, c = −1.765× 10 −5 , d = 9.843 × 10 −8 , e = 1.327 × 10 −8 and f = −297.Table 5 shows that for the ASTER DEM it was possible to reduce the error by almost two thirds (from 42.1 to 14.4), although there was still an off-set of 3.28.The standard deviation (STDEV) and average difference are computed over all the rocky areas higher than 4400 m (i.e.533 689 pixels).Errors calculated on exposed rock are not necessarily fully representative of the potential systematic errors on the glaciated terrain.
While this study was carried out in 2004-2005, a parallel study was ongoing using a similar approach on Coropuna, of which we were informed later (Racoviteanu et al., 2007).They built a DEM based on a topographic map of 1955, but at a more precise scale of 1:50 000.They used DEM from ASTER 2001 and SRTM 2000 datasets.Although the initial approach is similar (with exception that we used an additional DEM from 1997 and the ASTER 2002), there are several significant differences in the results of the two studies: Racoviteanu et al. succeeded to bring the standard deviation on rocky area to 9.5 m for SRTM (DEM 2000) as compared to 13.2 m in our study.This might be explained by the use of a more precise map (1:50 000 as compared to our 1:100 000), but could also be due to a smaller sample of reference points to compute the standard deviation (61 points in their study as compared with 533 689 in our present study for rocky areas).
However, we succeeded in decreasing the ASTER standard deviation to 14.4 m (in contrast to 26 m in their study), by using an additional correction based on a quadratic regression that corrected the ASTER 2002 DEM for bias induced by elevation and the cosine of aspect.The difference observed between DEM 2002t2 and DEM 1955 is −9.4 m (see Table 7), a result close to the results found between DEM 2000 and DEM 1955 (−8.75 m).Racoviteanu et al. (2007) found an elevation difference between 2001 and 1955 of + 28.5 m, which is extremely unlikely and not supported by both our and their GPS measurements.This is acknowledged in their article: "Comparison of GPS points with corresponding ASTER elevations on glaciated areas shows that the ASTER DEM is too high on glaciated terrain, with a RM-SEz error of 98.3 m with respect to GPS points".This highlights the importance for further correcting the DEMs with elevation and cosine of aspect.It also highlights that obtaining more precise maps should not be underestimated.

Change in ice cover
The identification of ice cover using images from passive satellite sensors is very straightforward in this location.As long as images are selected outside the snowy season (i.e.not between December and April), the ice detection is facilitated by the high contrast between the dark volcanic rock and the ice.The summits have a smooth shape, thus without much shadows (see Fig. 5).The images were projected in UTM 18S (datum WGS84), they were georeferenced, classified.
The classes corresponding to ice cover were transformed into vectors and the area computed using GIS.
The map in Fig. 2 shows the ice cover changes for the five dates.This first analysis revealed that the Coropuna Glacier shrunk steadily from 122.7 Km 2 in 1955 to 48.1 km 2 in 2008, i.e. losing more than 60% of its surface in 53 years (Table 6).
By plotting glacier area through time (Fig. 6) a clear declining trend appears.However with these observations it is not possible to predict whereas this will follow an accelerating trend (A) corresponding to smaller volume of ice having less inertia, thus shrinking faster.A linear trend (B), or a decelerating trend (C), where the shrinking will be slower when affecting higher altitudes.While scenarios A and B do not make much difference (total decline around 2040), in the  scenario C, a small glaciated area would be maintained at the higher altitudes and slowly decline.The three scenarios have a very good fit with R 2 of 0.996, 0.989, 0.987 for A, B and C respectively.A better understanding of the rate of volume loss and remaining ice volume, might provide additional insight toward which scenario is most likely.For example if we see a decline below a certain altitude but a steady (or increasing) volume, this might provide indication that scenario C is more likely than scenarios A and B. Conversely if the decline is all over the glacier scenarios A or B would be more likely.

Ice volume loss
Except the DEM 1997 (of which 24% did not contain data especially with respect to glaciated area), both SRTM 2000 and ASTER 2002 provide similar elevation differences.According to these results the average elevation changes were −8.75 ± 13.2 m (DEM 2000) and −9.4m ± 14.4 m (DEM 2002) at 95% confidence interval (see Table 7 for values and Table 5 for the margin of error).This corresponds to a yearly average loss of about −0.2m ± 0.3 m per year (at 95% confidence interval).Once extrapolated to the volume, the loss of ice between 2002 and 1955 is estimated to be 1.2 km 3 .This corresponds to an ice volume decrease of 18% as compared with 1955.These figures should be taken with caution given the low accuracy of the DEM, as illustrated by the measured difference being lower than the margin of error.
Unfortunately the DEM 1997 includes numerous pixels without data, especially on glaciated areas, making it less reliable.Figure 7 shows the significant improvements (especially on ASTER DEM) that can be achieved after correction of the DEMs.The margin of errors (see Table 5) is still higher than the measured differences (see Table 7), but this error follows a normal distribution with a clear higher amount of pixels depicting a decrease in elevation change (see Fig. 7).The much smaller differences observed in 1997 as compared with 2000 and 2002 are more likely due to data quality rather than linked with a 6 m decrease in ice thickness between 1997 and 2002.
Figure 8 shows the differences between the DEMs based on the topographic map from 1955.The 1997 DEM includes 24% of no data (in black) due to shadow of relief.The loss of thickness is represented in orange.The losses and the gains are located in the same areas, which is a good indication of the method's robustness, as these results come from three different DEMs generated from three different types of sensors.The differences in the south (at the bottom of the glacier) might be explained by seasonal changes.The SRTM mission (DEM 2000) was in February at a time with high snow cover, while the DEM 1997 was acquired in October, at the end of the dry season.DEM 2002 was acquired in May.

Analysing the transects
Figure 9 shows a portion (about 50%) of transect 2. This part is interesting as it shows that in some areas, the bedrock (in orange) is too deep to be recorded (see in the Crater).

Modelling remaining ice: results
In order to extrapolate the estimation of the ice volume to the rest of the map, the hypothesis was made that the depth of ice was dependent on the altitude, the aspect (slope's orientation) and the slopes.The assumption was made that ice thickness would depend on these three components.The amount of precipitation should be driven by altitude, slope and orientation.Orientation should also play a role due to predominant wind direction as well as different solar exposure (also probably less prominent in the tropics).Finally, snow accumulation should also be driven by slopes, based on the hypothesis that on steep slopes the ice should be thinner as the glacier is moving faster, whereas on gentle slopes, the ice accumulates as the glacier slows down.
To test these hypotheses, a statistical multiple regressions analysis was made using the recorded depth in relation to altitude, slope and orientation.
It was necessary to differentiate six cases (see Table 8): the top of volcanoes with altitudes higher than 6360 m; areas with altitudes between 6100 and 6360 m; and altitudes between 5980 and 6100 m.The next categories used three differentiations of slope orientation.These were needed for altitudes lower than 5980 m with the following aspects: higher than 270 • , between 91 • and 270 • and smaller than 91 • .Table 8 describes the variables (altitude, slope and orientation) and weight (in bold) used in the model according to the different thresholds.The quality of the models was assessed by looking at p-value1 (all smaller than 10 −10 ) and Pearson coefficient (between 0.80 and 0.94 for altitudes higher than 5980 m).The model becomes less accurate for lower elevations, which is reflected by a lower correlation (between 0.64 and 0.77 except one at 0.93).
Equations from Table 8 suggest that except for case 4 where orientation of slopes seems to play a role, in all the other cases, the depth of ice can be explained by altitude and slopes only.At the summit (altitude > 6360 m) ice depth is only linked with altitude.This is not surprising given that the smooth round summits of Coropuna are mostly flat (hence limited slopes and orientation).The map in Fig. 10 shows the result of this model once extrapolated to the entire glacier.
The modelled thickness ranges between 20 and 200 m, with an average thickness estimated at 80.8 m ±16.5 m (at 95% confidence interval).This gives an expected remaining ice volume of 4.62 km 3 (± 20.3%).The margin or error is  fairly important at ± 0.94 km 3 ; however, this result only provides a rough estimate.
The multiple regression analysis confirmed that altitude, slopes and orientation are factors influencing the ice thickness.This illustration includes transect 2, 3 and 4 (see Fig. 2 for their locations), thus offering the longest stretch across the measurements.The modelled depth fits well with the recorded profiles.The Pearson correlation between ice depths measured and modelled is 0.87 with a standard deviation of 16.5 (at 95% confidence interval).Figure 11 provides the fit between the observed and modelled depths.Still, for obvious reason of access, we were not able to take measurements on steep slopes with the GPR.This lack of samples might have an effect on the model.The maximum ice thickness on steep slopes, especially below the west summit, might be a limitation of our model to capture these physical conditions.The DEM 2000 shows an increase in elevation in higher than 6160 (areas where the green line above the blue line in Fig. 11).Elevation loss can be seen on lower summits (area where blue line is below green line in Fig. 11).

Estimation of the evolution in ice cover
This technique based on passive sensors provides simple and clear results.This is still the most effective method for quickly identifying a trend; however without remaining ice thickness and estimation of volume loss rate, it is difficult to www.the-cryosphere.net/4/313/2010/The Cryosphere, 4, 313-323, 2010 say whereas the projected decline will be abrupt or slow.It is possible that with a new estimation in year 2015, the trend might be better identified.

Measuring ice thickness
The hypothesis made on the statistical link between altitude, slope and orientation proved to be successful, except for orientation, which plays a limited role (mostly non-significant except in one case).It was expected that orientation would play a bigger role due to direction of the predominant wind, precipitation coming mostly from one direction.The limited role played by orientation is probably due to the study area's location within the tropics where the sun is less predominant on a specific slope (such as south in northern latitudes or north in southern latitudes).For glaciers at higher (Northern Hemisphere) or lower (Southern Hemisphere) latitudes, slope orientation might play a bigger role and should not be disregarded in the model.

Limitations on measuring the difference in ice volume
In this analysis, the differences recorded on the rocky areas were significant and raised concerns on the ability to use DEMs.In any case, they could not be used without the adjustments.The corrections applied on easting and northing as well as on elevation and aspects (for DEM 2002) significantly improved the quality of the DEMs.The Radar DEM 1997 included numerous no-data.Another limitation lies in the fact that digital DEMs are fairly recent technologies, and the date of the oldest archive found for the Coropuna was 1997.
The methodology for identifying the ice volume loss using DEM is straightforward and less challenging compared to the estimation of remaining ice.The difficulty comes from the lack of precision of DEMs on rough terrain, and attempts to correct the distortions can be time-consuming.

Conclusions and perspectives
The simple computation of area decrease based on passive satellite sensors, show a rapid decline of Glacier Coropuna (−1.4 km 2 per year, −60% since 1955) and if this continues, it may be gone before 2045.This is why it was important to try to produce more complex evaluations in terms of volume loss and remaining ice.
The methodology chosen for ice thickness and ice volume loss estimation proved to be effective.The corrections applied on the DEMs through multiple regression models based on easting, northing, aspect and elevation reduced the uncertainties, although the margin of errors is still high.
The vertical accuracy for the differences computed between the DEM2000t, the DEM2002t and the DEM 1955 were estimated at ± 13.2 and 14.4 m respectively, for elevation changes of −8.75 and −9.4 m on the ice (i.e. an average of 0.19 to 0.2 m of decrease per year, ± 0.3 m).Errors calculated on exposed rock are not necessarily fully representative of the potential systematic errors on the glaciated terrain.Although this provides an estimate of ice loss trends, the difference in ice thickness is smaller than the margin or error, thus affecting the confidence in the measurements.This can only be improved by either satellite sensors with better precision or by pursuing research in finding methods for correcting the DEMs.
For the ice thickness, the methodology could be improved, notably by choosing a lighter GPR.This would have eased the data collection, for instance by using skis to cover a much bigger area.Using the profiles from the ground study and a statistical extrapolation (modelling), it was possible to estimate the ice thickness (in average 80.8 m ±16.5 m), which gives an estimated remaining volume of 4.62 km 3 ±0.94km 3 (i.e.3.2 million tonnes of water).
These results were presented to GTZ and COPASA in Arequipa in December 2005.It helped to raise awareness on the issue of shrinking glaciers.In 2006, COPASA obtained support from GTZ.We then presented our results to United Nations Development Programme (UNDP) and they agreed to support COPASA through their Global Risk Information Programme (GRIP).Later on, the Inter-American Development Bank (IADB) also joined.These institutions help to introduce new policies for climate change adaptation at both local government and community levels: between 2006 and 2009, the following actions were carried out: A "Changing climate scenario" was developed for the Arequipa region; the socio-economic consequences of climate change were assessed.This led to a climate change adaptation strategy which was included (and implemented) in the Development Plans of six districts of Arequipa State; two urban and rural land use plans were developed in the Viraco and Machahuay districts; guidelines were developed with the Ministry of Agriculture for the incorporation of climate change adaptation in agricultural procedures; several thematic networks of students and teachers have been created and are working on climate change topic; the issue has been brought to the attention of regional institutions, which then produced a regional strategy for climate change adaptation.Three university theses have studied climate change at the regional level; four educational brochures were developed and their use approved in primary and secondary schools; a board game was developed on Climate Change Adaptation to help children to learn while playing.Five mini reservoirs and 15 warehouses for forage were built in this area.
The recent 2008 ASTER image shows that the glacier area continues to shrink, however, the local authorities have now integrated this threat into their development plan.The threat on water supply might be increasing, but efforts are made to reduce the vulnerability of the local population.
1) Where I = Ice thickness [m], T = Time [ns], C = Speed of propagation through ice of the signal (0.16 [m/ns]) For example a time lag of 2000 ns = 160 m of ice thickness.
shows that they are not in the same plan, with DEM 1997 having the least distortions and DEM 2002 (ASTER) showing the worst distortion.If the DEMs are not in the same plan, a change in latitude and longitude can significantly influence the difference in elevation

Fig. 4 .
Fig. 4. Difference on rocky areas versus elevation and aspect for ASTER 2002.

Table 4 .
Weights used to place the DEMs in the same plan as DEM 1955.

Fig. 7 .
Fig. 7. Elevation differences on rock and on glaciated areas.

DEM97
40 m * * Once the off-set of 3.28 is removed, otherwised the average difference on ice is −6.11.

Fig. 9 .
Fig. 9. Example of interpreted profile of transect 2 from the GPR on the slope to the summit.

Table 1 .
Data source for monitoring glaciated area.

Table 2 .
Sources and general information for the digital elevation models.
Fig. 2. Georeferenced radar profiles and evolution of glaciated area

Table 5 .
Assessing the DEMs accuracy before and after corrections (on rocky areas).

Table 6 .
Evolution of ice cover.

Table 7 .
Altitude changes on glaciated areas on original and corrected images.

Table 8 .
Quality of regressions used for modelling ice thickness.