Assessing high altitude glacier volume change and remaining thickness using cost-efficient scientific techniques : the case of Nevado Coropuna ( Peru )

Higher temperature and change in precipitation patterns have induced an acute decrease in Andean glaciers, thus leading to an additional stress on water supply. To adapt to climate changes, local governments need information on the rate of glacier volume losses and on current ice thickness. We show how volume changes can be estimated in remote areas using readily available low-cost digital elevation models derived from both topographic maps and satellite images. These have been used for estimating the volume changes over the Coropuna glacier (Peru) from 1955 to 2002. Ice thickness was measured in 2004 using a georadar coupled with 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 regressions techniques. Computers were modified to resists high altitude (6500 m) temperatures and low pressure conditions. The results delineated a significant glacier volume loss and provided an estimate of the remaining ice. It provided the scientific evidence needed by local Peruvian NGO, COPASA, and the German Cooperation Program in order to alert local governments and communities and for enforcing new climate change adaptation policies.


General context
Changes in glaciers and ice caps are good indicators of evidence of climate change (Zemp et. al., 2008).A majority of the world glaciers have undergone a reduction in their mass at an accelerating rate (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 (Stern, 2007, p.56).
In Peru, the demographic 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.Melting glaciers increase flood risk during the wet season and reduce dry-season water supplies.This is of concern in the Indian sub-continent and parts of China as well as in the Andes (Stern 2007, p.56).In the latter 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 600 km 2 in 1996 (Silverio, Jaquet, 2005) in Cordilera Blanca (Peru).

Assessing change of ice volume in Nevado Coropuna (6500 m, Peru)
The present study on Coropuna Glacier was made at the request of the Cooperation Peruano Allemana Para la Seguridad Alementicia (COPASA, Special Project of Arequipa Regional Government in collaboration with the German cooperation program (GTZ).They needed scientific evidences of glacier retreat in order to introduce climate change adaptation policies on water supply along with local government and communities.
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.
A first analysis based on archive satellite images (Silverio 2005) revealed that the Coropuna Glacier lost 54% of its surface between 1955 and 2003.It also showed that the level of precipitation varies significantly during El Niño events.This paper concentrates on the second part of the study, which was to measure the change in volume of the glacier using different Digital Elevation Models (DEM).In order to measure the current (2004) ice thickness, a field mission was carried out using a georadar coupled with Ground Positioning System.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 and statistical multiple regressions techniques.
Given the limited financial resources of the local governments, and development organisations, simple and low-costs techniques to follow glaciers volume dynamic and remaining ice volume are needed.The purpose of this paper is to show how with limited budget and using locally available hardware, scientific evidences of glacier volume variation and ice thickness can be obtained.

Study area
The Coropuna Glacier is the third highest summit in Peru, culminating at 6,446m.It is located at 15.546 S, 72.660 W, about 155 km northwest of the city of Arequipa in Peru.According to COPASA staff, 8,000 persons depend on the Coropuna Glacier for their water supply and it is estimated that 30,000 people depend indirectly on the glacier for their livelihood.
The study area covered 577 km 2 around the Coropuna summit.This area includes the entire glacier surface as observed in 1955.A base camp was set up at 5664 m, and ground measures were collected at altitudes ranging between 5,870 and 6,446 m (see Figure 1).

Digital Elevation Models
In order to estimate the ice volume loss between 1955 and today, 8 Digital Elevation Models (DEM) from different years and periods were considered.The first one was generated with manually digitalised isolines from the topographic map of 1955.The one from 1997 was acquired by GTZ/COPASA from SARMAP, based on a pair of ERS radar satellite images; the DEM from 2000a was based on radar measures from the space shuttle.All the others were based on ASTER satellite data (see Table 1).This choice was made based on the quality of the dataset, but also to ensure an adequate time span between the datasets.It should be noted that the DEM 2003 would have offered a better time span, but since only half of the glacier was covered, it was finally rejected.

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 DEM.The 14 day-mission was undertaken between 13 and 26 August 2004.The team was composed of two scientists and 11 support staff (guides, porters,…).
The scientific instruments used were selected based on local availability.The Ground Penetrating Radar (GPR) Ramac X3M, included a 100 MHz shielded antenna, batteries and electronic device.Much lighter GPR exists; however, we used the only GPR available in Arequipa.To avoid carrying the 35 kg antenna every day, it was protected by a waterproof bag and buried under the snow.Three Global Positioning System receptors (GPS) were used and two regular office laptops.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 measures were recorded on USB cards.

Estimation of ice thickness
Measuring the ice volume was achieved with the GPR.The device was set to emit at 438 MHz.The signal travels through ice at 0.16 m/ns.This setting was used according to other studies performed in similar conditions (Gruber, Ludwig and Moore, 1996).The depth of ice can be measured by recording the time lag between the emission and the reception of the signal (see Equation 1) Technical settings are specified in Table 2.This allows the detection of the bedrock.Figure 2 illustrates the set up used during the collection of the data.Each recorded depth was coupled with geographical coordinates obtained using a GPS so that the profiles could be geo-referenced.Given time and access constraints, it was not possible to

Record parameters Settings
Sampling frequency 438 MHz

Number of samples 900
Number of stacks 16

Time window 2055 ns
Trace interval 2 seconds Antenna separation 0.5 meter achieve a comprehensive recording of the glacier surface.Instead, the mission proceeded along transects (shown on Figure 1).The selected targets were chosen to provide samples including different altitudes, slopes, and aspects.The hypothesis was that these three variables would explain most of the ice thickness.Using multiple regressions analysis, depth was modelled in areas where no measures were taken.
During the mission, 10.6 km of transects were realized using the GPR.The map in Figure 1 shows the distribution of GPS points recorded along these transects.A GPR signal was recorded every two seconds, and a GPS location with hours, minutes and seconds approximately every 2 minutes.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 detect ice depth.Due to the computer configuration that limited record time window, the bedrock was sometimes too deep to be detected (typically in volcanoes craters, see Figure 3).However, it was possible to extrapolate by following ice stratifications.Whenever this was possible, such evaluations were processed.

Estimation of ice volume loss
Although passive satellite sensors, such as Landsat TM, provides estimate of the area covered by ice (Silverio, Jaquet, 2005).It does not provide information on the loss of ice thickness.
The idea is to measure the difference of ice volume using DEM time series.

Adjustment and corrections
The vertical accuracy of available digital elevation models remains very low.The free DEMs provided by ASTER have a nominal vertical accuracy of 30 m according to the USGS, who generates them.However, studies made using these DEMs, reveal significant errors when applied to rugged terrain and steep slopes (Kääb et al. 2002).
In order to compare different DEMs, several operations were needed.Firstly, all the DEM were re-sampled to 30 m to compensate for different spatial resolutions.They were then reprojected in Universal Transverse Mercator (projection 18 south, datum WGS84).Finally, they were geo-referenced so that they could be overlaid.This was performed using ground control points, such as summits located outside the glacier area (on bare rocks).
An adjustment for vertical accuracy was required given that the DEMs had different altitudes for over the same points on rocky terrain (-46m to +53m as compared with the map 1955).
This was achieved by using a simple vertical off-set, using rocky points as references.This is an approximation, given that radar images are behaving in different ways on different landcover (Liu, 2006;Rabus et al, 2003).Several other adjustment methods were tested: ranging from simple adjustment to modelling of errors using linear regression or multiple regressions.A reference area was chosen on surfaces that were not covered by ice, snow or vegetation, i.e. rocky, bare ground on all DEMs.This corresponds to the area outside the glacier area of 1955.The hypothesis was that observed difference on a reference surface could be modelled using slopes, altitude and orientations.These factors proved to be influencing the DEM errors in other studies (Gorokhovich and Voustianiouk, 2006).Although statistical methods provided good results on the reference surface, it may have introduced some artefacts in a dynamic glaciated area and finally, only a simple off-set method was used.The issue of distortion was especially significant on ASTER images, where extensive correction would be needed (altitude are tilted with X and Y).
While this study was carried out in 2004 -2005, a parallel study was ongoing using similar approach on Coropuna, of which we were informed later (Racoviteanu et. al, 2007).They also used the 1955 map as reference for the DEM as well as ASTER and SRTM datasets.Their Once the correction was applied, the average altitude on bare rocks was equal between all DEMs.Over the ice the sum of all altitudes was then computed (all pixel values) for all the DEM and the sum of all altitudes for 1955 was subtracted from this value.The result is divided by the number of pixel to get the average difference between the DEM and the status in 1955 (see Equation 3).Using the sum of all the altitudes (sum of all pixel values) over the whole reference area (rocky, bare ground) instead of a pixel by pixel subtraction in a raster GIS allowed avoiding the errors introduced by the geo-reference.Indeed, in mountainous areas, a 60 m horizontal difference in geo-location can lead to several hundred metres of vertical difference.Summing up all the individual pixel by pixel differences, will introduce errors from the resampling method.Since DEM were originally at different spatial resolutions, it is suggested that matching individual pixel is not appropriate.

Analysing the transects
Figure 3 shows a portion (about 50%) of analysed 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).In this area, the bedrock is extrapolated following lines of ice stratifications.

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 slope orientation (aspect) and the slopes.The assumption was made that the accumulation of ice would be somehow affected and depending from these three components.Amount of precipitation should be driven by altitude, slopes and orientation.Orientation should also play a role due to predominant wind direction as well as different solar exposition (also probably less prominent in the tropics).Finally snow accumulation should also be driven by slopes, making 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 profile recorded depth in relation with altitudes, slopes and orientations.
It was necessary to differentiate six cases (see Table 3): the top of volcanoes with altitudes higher than 6360, then areas with altitudes comprised between 6100 and 6360, then for altitudes between 5980 and 6100.The next categories used three differentiations of slope orientation.These were needed for altitudes lower than 5980 (aspect higher than 270°, between 91° and 270° and smaller than 91°).Table 3 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 were 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).The model becomes less accurate for lower elevations, this being reflected by a lower correlation (between 0.64 and 0.77 except one at 0.93).Equations from Table 3 suggest that except for case 4 where orientation of slopes seems to play a role, for all the other cases, depth of ice can be explained by altitude and slopes only.
On the summit (altitude > 6360) ice depth is only linked with altitude.This is not surprising given that the smooth round summits of Coropuna have limited aspect variation and are mostly flat (hence limited slopes and orientation).The map in Figure 4 shows the result of this model once extrapolated to the entire glacier.The thickness ranges between 20-200m, with an average thickness estimated (based on the model) between 91.6 and 95m.If extrapolated to the whole area, the volume is expected to be between 5.25 and 5.46 km 3 .
The multiple regression analysis confirmed that altitude, slopes and orientation are factors influencing the ice thickness.Figure 5 shows the fit between the model and the profile 2. This illustration includes transect 2, 3 and 4 (see Figure 1 for their locations), thus offering the longest stretch across the measures.The modelled depth was found to follow the recorded profile with good fit.Still, for obvious reason of access, we were not able to take measures on steep slopes with the GPR.This lack of samples in steep slopes 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 apprehend these physical conditions.
The results presented in these profiles are consistent with what was expected, except maybe in transect 4, where the ice was thicker than expected in North-East direction (between 6000 and 6100 m).However, at higher altitudes and with other orientations, the results provided a good fit between the model and the measures.Figure 6 shows the location of the different transects across the model, while Figure 7 depicts the variation of ice thickness (in blue) as modelled across the five different profiles.The findings show that the ice is thicker at high altitude and on gentle slope.

Ice volume loss
The "simple adjusted method" as used in Equations 2 and 3, was applied.Table 3 provides the average elevation on the reference area (RA) and on the glacier (GL).An average thickness difference of 16.7 m was found by subtracting average altitude on the glacier in 1955 with the same areas in 1997 (see Table 3).A difference of 10.5 m was found using the same process applied between 1955 and 2002.Although the Coropuna region received much precipitation during the El Niño 1997/98, the smaller differences observed in 2002 as compared to 1997 is more likely to be due to data quality rather than linked with a 6 m increase in ice thickness between 1997 and 2002.The 1997 radar DEM is more precise on the rocky area, although it includes a significant amount of "no data" findings.The margin of error being important, it is difficult to guess which of the DEM difference is closer to the reality.According to these results the average yearly losses could vary between 0.22 and 0.4 m per year.Once extrapolated to the volume, the loss of ice between 2002 and 1955 is comprised between 1.2 and 2.1 km 3 .This corresponds to a decrease varying between 17.9% and 28.8% as compared with 1955.These figures should be taken with caution as the DEM has low accuracy.

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-significative except in one case).It was expected that orientation would play a bigger role due to direction of predominant wind, precipitation coming mostly from one direction (northeast),...The limited role played by orientation is probably due to the location within the tropics where the sun is mostly vertical.For glacier at higher (north hemisphere) or lower (south hemisphere) latitudes, slopes orientation should play a bigger role and should not be disregarded in the model.
The accuracy of the model was better at high altitude than at lower altitude (< 5980).This is believed to be due to two factors.Firstly there were more records made at high altitude than at lower altitude (due to difficulty of access on the edge of the glacier), secondly, it is believed that at lower altitude, site effects played a more significant role hence were more complex to model.This could be due to lower temperature inertia (thinner ice) and increased solar heat radiation from the surrounding rocks.

Limitations on measuring the difference in ice volume
In this analysis, errors on the rock were significant and raised concerns on the ability to use DEMs (especially ASTER).In any cases, they could not be used without the adjustment.
Radar images proved to be more precise than the free ASTER DEM, however Radar sensor are fairly recent technologies and the date of the oldest archive found 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 DEM on rough terrain, where attempt to correct the distortions appeared to be quite time consuming.It is expected, however, that such technologies will evolves as numerous new radar sensors are now available, including some with finer spatial resolution.

Conclusions and perspectives
The methods chosen for ice thickness and ice volume loss estimation proved to be efficient, although the choice of a lighter GPR 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 between 91.6 and 95m), which gives an estimated remaining volume comprise between 5.3 and 5.5 km 3 (i.e.3.7 and 3.8 millions tons of water).
The study using the DEM accounted for most of the imprecision.DEM elevation accuracy is estimated to be +/-30m and statistical corrections were needed to make the variations DEM versions comparable.Nevertheless, it shows that the Coropuna glacier is diminishing both in area and volume.The average thickness loss is estimated between 0.22 and 0.40 meters per year.
It is not possible with the current number of DEMs to say if the loss is linear or if it will accelerate.To assess this, new measures would need to be taken, every 4-5 years.It is recommended to use radar satellite images.The more expensive radar DEM proved to be more accurate than ASTER passive sensors, although it included numerous no data pixels.If DEM are taken every five years or so, within 15 years it should be possible to say if the ice reduction is linear or is accelerating.
In order to avoid the risk of including snow cover in the estimates, the best period for data acquisition is during the months of September -October.
The results from this study were presented to GTZ and UNDP.Although the ice volume loss was strongly suspected, having factual numbers to present to these development agencies, helped to convince them to continue supporting research and irrigation projects.A follow-up study for evaluating the impacts on water supply might consist in assessing where the most important water sources are located and how they might be impacted by the glacier retreat.A precise delimitation of water catchments, measure of river flows at different places are tools that would allow a refined planning for future irrigation of crops.
of propagation through ice of the signal (0.16 [m/ns]) For example a time lag of 2000 ns = 160 m of ice thickness.
study used a pixel by pixel subtraction, which, given the low level of precision and rugged terrain, the authors of this paper felt that such an approach requires a level of precision that exceed what ASTER DEM can provide.The conclusions of Racoviteanu et.al's mentioned that "Spurious values on the glacier surface in the ASTER DEM affected the analysis".It potentially prevented them from quantifying the glacier changes based on the ASTER DEM.The same problem was encountered by the study described in their paper.An adjustment was necessary.It required computing the average altitude on reference areas covered by rocks on all DEMs.The altitude of reference was chosen using the DEM 1955.A first process was applied by adjusting the difference in altitude observed on bare rocks, by subtracting the difference with the DEM of reference (see Equation 2pixels on bare rock (area defined based on bare rock as in 1955) for the DEM "A" P BR 1955 Altitude value of pixels on bare rock for the DEM 1955 N BR 1955 = Number of pixel in bare rock as of zone corresponding to 1955 Average difference of ice thickness between "DEM A" and "DEM 1955" P ice A= Altitude value of pixels (area define based on ice in 1955) for the "DEM A" P ice 1955= Altitude value of pixels (area define based on ice in 1955) for the DEM 1955

Figure 3
Figure 3 Example of interpreted profile of transect 2 from the Georadar on the slope to the summit.

Figure 4
Figure 4 Estimation of ice thickness (model)

Figure 5 Figure 6 Figure 7
Figure 5 Comparison between ice thickness measured and modelled over transect 2, 3 and 5

Figure 8
Figure 8 Estimation of ice loss between 1955 and 2000 and between 1955 and 1997 Figure 8, shows the differences between DEM based on the topographic map 1955 and the DEM from SRTM 2000 (map on top) and the radar images 1997 (map below).The DEM includes 24% of no data (in black) due to shadow of relief.The loss of thickness across timeis represented in orange.The same computation using the ASTER DEM presented a smaller average difference, although with higher variations and seems less reliable as compared to radar DEM 1997.The SRTM is provided here to show the seasonal variation.The SRTM mission was in February hence at time with high snow cover, while the DEM 1997 was acquired in October, at the end of the dry season.This explains part of the difference in thickness.