Zhen Wu,ShiYin Liu,ShiQiang Zhang,HongLang Xiao

1.Key Laboratory of the Ecohydrology of Inland River Basins/Heihe Key Laboratory of Ecohydrology and Integrated Rivers,Cold and Arid Regions Environmental and Engineering Research Institute,Chinese Academy of Sciences,Lanzhou,Gansu 730000,China

2.State Key Laboratory of Cryospheric Science,Cold and Arid Regions Environmental and Engineering Research Institute,Chinese Academy of Sciences,Lanzhou,Gansu 730000,China

Englacial hydrological characteristics of a typical continental-type glacier in China as detected by ground-penetrating radar

Zhen Wu1*,ShiYin Liu2,ShiQiang Zhang1,HongLang Xiao1

1.Key Laboratory of the Ecohydrology of Inland River Basins/Heihe Key Laboratory of Ecohydrology and Integrated Rivers,Cold and Arid Regions Environmental and Engineering Research Institute,Chinese Academy of Sciences,Lanzhou,Gansu 730000,China

2.State Key Laboratory of Cryospheric Science,Cold and Arid Regions Environmental and Engineering Research Institute,Chinese Academy of Sciences,Lanzhou,Gansu 730000,China

The englacial structures and ice thickness of the Laohugou No.12 (L12) Glacier in the Qilian Mountains,China,were retrieved from ground-penetrating radar (GPR) profile data acquired in August of 2007.Here the interpretation of a typical GPR image is validated using two-dimensional,Finite-Difference Time-Domain (FDTD) numerical modeling.Data analyses revealed many englacial characteristics,such as temperate ice,crevasses,and cavities at the position of convergence between the eastern and western glacial branches of L12,and at an altitude between 4,600 and 4,750 m a.s.l.on the east branch.Combining ice thickness,englacial structures,subglacial topography,and surface flow velocities of this glacier,we analyzed the reasons for the distribution of temperate ice.The results show that greater englacial water content is associated with englacial crevassing and surface moulins,which allow water to be channeled to the temperate ice aquifer beneath the surface cold ice layer.Analysis of air temperature data shows that as more meltwater imports into the ice body,this has a great effect on water conservation and dynamics conditions.With climate warming,and under the influence of crevasses,subglacial structures,and ice thickness,ice thickness reduction on the L12 east branch is more rapid than that on the west branch.

ground-penetrating radar (GPR);Finite-Difference Time-Domain modeling;temperate ice

1.Introduction

The advance and retreat of mountain glaciers are sensitive functions of climatic change (Su and Shi,2002;Barry,2006).The Earth’s climate has always been characterized by natural variations.However,mean annual air temperatures rose dramatically in the 20th century (Shi and Liu,2000;Alleyet al.,2007;Molnia,2007).This trend intensified at the end of the last century,and glacial ice coverage greatly diminished in the Alps and in central Asia (Paulet al.,2004;Khromovaet al.,2006).Glacier retreat has not been homogeneous;it has depended strongly on glacier size,location,and climate regime.For example,in the 20th century,glaciers in the more continental Pamirs retreated less than those in the more humid parts of the Tianshan Mountains (Liuet al.,1999;Hagg and Braun,2005).

Similarly,for individual glaciers,retreat has also not been homogeneous,but has depended strongly on air temperature,englacial structures,hydrological characteristics,and climate regime (Irvine-Fynnet al.,2005).With greater glacial ablation,englacial crevasses and cavities with water have increased (Murray and Porter,2001;Jørgensen and Andreasen,2007).The rheology of ice (ice softness) varies with water content;a wet ice bed is more likely to allowsliding and bed deformation (Iken and Bindschadler,1986).The variation of dynamic characteristics has a particularly great impact on mountain glaciers (Serranoet al.,2006;Golledge,2007;Hooke and Fastook,2007).Staging saltation is linked with temporary changes of hydrological characteristics (Bowersox,2002;Herzfeld and Mayer,2003;Thorsteinssonet al.,2003).The detection of hydrological characteristics is usually obtained from ground-penetrating radar (GPR).

Many researchers have used GPR on mountain glaciers and polar ice caps to survey englacial structures and ice thickness (Pettersson,2004;Harperet al.,2005;van de Veenet al.,2007).For instance,Taurisano (2006) and Pattynet al.(2009) investigated the distribution,size,and causation of crevasses in Antarctica;others surveyed the distribution of ice thickness (Thorsteinssonet al.,2003;Binghamet al.,2005;Tamuraet al.,2006;Welchet al.,2009;Shugaret al.,2010;Bartholomauset al.,2011;Haaset al.,2011),or simulated and analyzed the evolution of drainage systems and showed the effects of glacial ablation on englacial structures (Thorsteinssonet al.,2003;Binghamet al.,2005;Bartholomauset al.,2011).

Radio echo soundings and GPR surveys of Arctic glaciers have identified a spatially extensive,internal reflecting horizon in several glaciers and ice sheets.This feature is commonly interpreted as either the cold ice or warm ice glacial boundary (Wohllebenet al.,2009).Warm ice is defined as ice enriched with water;the existence of warm ice has been verified by comparison of GPR data and ice core samples (Fritzscheet al.,2005;Davis and Thompson,2006).The warm ice zone has been found in many continental glaciers and polar ice caps in recent years (Mooreet al.,2009;Humbert,2010).

Although the discovery of the warm ice zone has been common,there has been less study on its causes and its indication of climate change.In this study,we combine GPR and other data,and analyze the cause of warm ice from several aspects of subglacial topography,ice thickness,glacial velocities,distribution of crevasses,and glacial morphology.These results provide a reference for future research into the mechanisms of glacial response to climate.

2.Glaciological setting and sampling methods

2.1.Site description and equipment

The L12 Glacier (Figure 1) is the longest valley glacier located in the Qilian Mountains (39°25′–39°30′N,96°30′–96°35′E);it is composed of two separate branches,the east and the west (Sun and Xie,1981;Duet al.,2008).This glacier is a typical extreme continental type (Wanget al.,1981).The L12 terminus is at 4,260 m a.s.l.,and its highest altitude is 5,483 m.There are numerous seracs from the terminus to an altitude of 4,500 m,with a wide firn accumulation region.The slope of the entire glacier is gentle(3°–6°),and there are many supraglacial rivers and melting caves on the surface.

Figure 1 Location and schematic presentation of the Laohugou No.12 Glacier

We investigated the L12 on August 1 and 8,2008,using a pulse EKKO PRO GPR (Sensors &Software Inc.,Mississauga,Ontario,Canada).The investigation time was the key glacial ablation period,as our investigation could better reveal hydrological characteristics during that time.Profiles were obtained using a 50 MHz bistatic antenna.The transmitting-receiving antennae were arranged parallel to each other with 4-m separation,transverse to the profile direction.The antenna step distance was 1 m to ensure adequate horizontal resolution.

The GPR profiling was run in a common offset mode and traces were triggered artificially.

Data were collected over an approximately 18 km stretch,between 4,400 and 5,000 m a.s.l.Profiles were recorded at one trace per meter,and each trace was resolved using 64-bit samples,which ensured the horizontal resolution of the measurement section.A global positioning system (GPS) device from Beijing UniStrong Science &Technology Co.,Ltd.(Beijing,China) was used to position each profile.The standard error of the GPS devices was 0.01 m ±1 per million (horizontal error) and 0.02 m ±1 per million(vertical error) under the survey and signal reception conditions.However,surveys of surfaces in RTK differential mode typically produce horizontal and vertical survey errors of 0.10 to 0.30 m,respectively,for geodetic-quality GPSreceivers (Riveraet al.,2005).These errors are acceptable for GPR data processing and analysis.

Radio wave velocities were calculated from the travel time to a point diffractor (Macheretet al.,1993) as observed in the radar record.We could use only symmetric diffraction hyperbolae,which corresponded to either point or horizontal linear diffractors that were crossed by the radar profile(Rolstad and Norland,2009).Features such as water-filled channels were detected and measured using multiple profiles to accurately estimate radio wave velocity (Macheretet al.,1993).At the same time,compared with other results of the continental glacier (Binderet al.,2009;Berthier,2010;Shean and Marchant,2010),we obtained a radar wave velocity of 0.167 m/ns,this result is consistent with previous study.

2.2.Numerical modeling

To test our interpretation of the englacial structures,which had not been directly explored and whose signatures appeared ambiguous,we used a Finite-Difference Time-Domain (FDTD) model (Joseph and Taflove,1997;Bergmannet al.,1998;Reineixet al.,2001;Wang and McMechan,2002;Lampe,2003;Liu and Zeng,2007;Guan and Hu,2008) that solves Maxwell’s equations in two dimensions (Figure 2).

Figure 2 Finite-Difference Time-Domain (FDTD) view of the model space

The same modeling technique has been used by Brandtet al.(2009) and others.The model was run repeatedly,to account for different geometries and media.Input parameters were based on each radar image to be validated.

For numerical stability and reliability,only features larger than 5 times the grid cell size were used,which were the smallest sizes for features that can give a realistic response.This size is coarse,but sufficient to produce a realistic pattern.The grid size was determined to be larger than 10 times the wavelength (λ) of the highest significantly powerful frequency in the pulse spectrum.A rule of thumb for a Ricker wavelet source is:

wherefmis the highest frequency in the pulse spectrum andfcis the center frequency.The wavelength (grid cell size) is then determined as:

wherecis the speed of electromagnetic waves in vacuum andεris the relative permittivity.The central frequency was set to 50 MHz and a grid cell size of 0.125 m was used in all simulations.The relative permittivity is related to density as described above (Bongiovanniet al.,2008).These were run under Perfectly Matched Layer (PML) boundary conditions.It should be stressed that FDTD simulation focuses on form rather than scale.

3.Results

3.1.Square cavity

Figure 3 presents a kind of englacial structure with three different media (air mixed with water,no water,and saturated with water).Overall,these three simulations show a downward hyperbolic reflection,but there are some differences in detail.

For an air cavity mixed with water,there is a largest dielectric constant difference between the two media,bringing about an obvious reflex arc at the air-water interface.

Because electromagnetic wave velocities vary within different media,when these waves pass through the boundary between two media,a velocity that changes from a variable state to stability would indicate some disturbance.Several superimposed hyperbolic reflections in Figure 3b1are related to velocity variations.Another possible explanation is that the strong extending resonant oscillation depends not only on the diameter of the void,but also on its height,or dimensional ratios (Kofmanet al.,2006;Wuet al.,2009).

For an air cavity (Figure 3a2),electromagnetic wave intensity shows weak reflection at the upper and lower interfaces of air and ice (Figure 3a3),indicating that electromagnetic waves pass through air with less energy loss.

For a cavity saturated with water,the reflection image is similar to hyperbolae from point sources.The result shows that electromagnetic waves have poor resolution for water.

Figure 3c shows typical hyperbola reflection,superimposed by only two hyperbolae.Reflection intensity is also greater.The model provides a radar image (Figure 3b3) of the cavity that is very similar to the observed signature (Figure 3c).For a glacier without supraglacial debris,strong reflection is only caused by englacial water.

Figure 3 Modeling of englacial rivers with the same form and different media distribution.Modeled geometry (air mixed with water)(a1);modeled geometry (no water) (a2);modeled geometry (saturated with water) (a3);modeled radar image (b1,b2,b3);actual radar image (EZ segment in Figure 1,elevation between 4,613–4,657 m) (c).

3.2.Crevasses

Figure 4 presents three simulation results for vertical crevasses of different forms.Figure 3a1presents observed and modeled radar images of a crevasse with vertical walls and two small internal bridges.The modeled radar image shows clear superimposed hyperbolic reflections,which suggests that the steep incidence angle of the radar signal on the vertical walls causes more energy to be deflected and lost downward and sideways.In this study,we cannot interpret the reason for this geometrical effect(scale effect).

In Figure 4c,there are several superimposed hyperbolic reflections on the glacial surface.Comparison of the simulation results indicates that the reflections are caused by crevasses in the adjacent glacial surface.According to the variation distance described above,it can be concluded that crevasse widths gradually decrease from top to bottom.

3.3.Temperate ice

We modeled within-ice water distributions at different scales (Figure 5).The simulation result shows that when the diameter of the origin (water) in the model exceeds the radar resolution,the radar reflection images show normal point-source reflection characteristics.The linear reflection in Figure 5b1is the result of point-source reflection superposition.When the diameter of the origin (water) in the model is less than the radar resolution (Figure 5a2),the reflection image shows a significant scattering phenomenon (Figure 5b2).Figure 5a3shows a different scale than Figure 5a2.The diameter of the origin (water) in Figure 5a3is less than that in Figure 5a2,but the simulation results in Figure 5b3shows the same image as those in Figure 5b2.This indicates that when the scale of water is less than the radar resolution,water with the same distribution and difference scale usually produces the same reflected image.

This scattering described above is very similar to the measured image of Figure 5c.Because of scattering in the lower ice body and the scattering area thickness reaching 40 m,it is unlikely that the scattering phenomenon results from sediment or englacial debris;it is almost certainly related to characteristics of the ice body itself.The simulation shows that the scattering phenomenon is highly dependent on the englacial small-scale water distribution.Therefore,it can be confirmed that this phenomenon is caused by water content rather than sedimentary features.

Figure 5c shows a dense chaotic return zone.Generally,when electromagnetic waves encounter materials with particle size less than the radar resolution,the waves will not be reflected.According to Kotlyakov and Plewes’s theory: (1) Water is known to cause perturbations in the wave profile (Kotlyakov and Macheret,1987);and (2)during the ablation season,scattering of electromagnetic pulses can result from enhanced melting and/or refreezing,leading to englacial heterogeneities (Plewes and Hubbard,2001).It may be concluded that the dense chaotic return zone is caused by an ice body with water,similar to the findings of Björnssonet al.(1998),Petterssonet al.(2004),and Irvine-Fynnet al.(2006).Englacial water content has been associated with ice beneath surface crevassing and moulins,which allows water to be channeled to the temperate ice aquifer beneath the surface cold ice layer (Mooreet al.,2009).

4.Distributions of typical characteristics

Through several investigations and verifications,we describe the distributions of typical characteristics for the entire survey trace (Figure 6).Englacial cavities with water are distributed at lower altitudes on the east branch,and at altitudes below the intersection of the east and west branches.Temperate ice is widespread below 4,550 m a.s.l.and at lower altitudes on the east branch,collocated with the englacial cavities with water.We also investigated the distribution of surface rivers and lakes on the glacial surface.Supraglacial rivers are mainly found below 4,600 m a.s.l.Supraglacial crevasses are mainly located on both sides of transverse section and at 4,600–4,700 m a.s.l.on the east branch.This reflects the position of maximum shearing strength near both sides of the bedrock,but not at the glacial edge (Reid and Callender,1965).With the expansion of glacial surface rivers,surface water transfers heat to adjacent ice by crevasses,thus contributing to the warming of the ice body in these regions.This drainage promotes the development of subglacial water systems and temperate ice.

5.Discussion

5.1.The effects of surface velocity on temperate ice

Supraglacial crevasses and englacial cavities with air usually result from uneven pressure caused by surface flow velocity differences.We obtained surface velocity in 2009–2010 mass balance years (Table 1).Comparison of the position of the ablation stake in Figure 6 shows there is little difference in surface flow velocity between position G of the wide firn basin and position E of the narrow transverse section,and the thickness at position G is greater than that at position E.If positions G and E have the same ice flux,the surface flow velocity at position E would be obviously greater than at position G.Nevertheless,actual measurements show that the surface flow velocity difference between these two positions (E (7.3 m/d) and G(7.8 m/d)) is not obvious.The ice body near altitude 4,700 m a.s.l.is likely being compressed by the upper ice body.However,the three large crevasses on the east branch between 4,650 m a.s.l.and 4,780 m a.s.l.are also probably being compressed by the upper ice body,generating surface uplift.

Figure 4 Modeling of continuous square crevasses.Modeled geometry (continuous same square crevasses) (a1);modeled geometry(continuous difference square crevasses) (a2);modeled geometry (long square crevasses) (a3);modeled radar image (b1,b2,b3);actual radar image (WZ segment in Figure 1,elevation between 4,550–4,650 m) (c).

Figure 5 Modeling of temperate ice.Modeled geometry (water point is greater than the radar resolution) (a1);modeled geometry (water point is less than the radar resolution) (a2);modeled geometry (large-scale water point is less than that in Figure a2) (a3);modeled radar image (b1,b2,b3);actual radar image (HZ segment in Figure 1,elevation between 4,420–4,500 m) (c)

Figure 6 The distribution of englacial and subglacial typical characteristics

Table 1 Surface flow velocities for the different ablation stakes in Figure 6

5.2.The effects of ice thickness and glacier bed on temperate ice

Ice thickness has an important effect on glacier surface velocity.The GPR results show that the average ice thickness on the east branch (181.1 m) is much greater than that on the west branch (130.8 m) (Figure 7a).There is little difference in surface slope between the branches.Therefore,the east branch should have greater ice flux,the result of which is likely to produce supraglacial and englacial crevasses.

As shown in Figure 7b,the topography of the glacier bed on the west branch may transfer the ice flow directly into a low-altitude ice body form accumulation zone,resulting in a more stable ice flow velocity.At an altitude of 5,000 m a.s.l.on the east branch,however,the ice flow initially concentrates near the center of the accumulation zone.Then the ice flow begins supplying the low-altitude ice body,recharging ice body in a short time.This supply on the east branch supports the production of both supraglacial and englacial crevasses.

5.3.The effects of meltwater on temperate ice

The crevasses and moulins allow the gradual import of water into the cold ice layer (Figure 8),increasing englacial water content under the influence of gravity (Petterssonet al.,2007).Meanwhile,the storage and transmission of englacial water have large effects on englacial structures and thermal characteristics (e.g.,meltwater may promote the development of new englacial waterways and broaden existing waterways,further affecting glacial melting rates),asthe rheology of ice (ice softness) is proportional to water content.As the englacial water content near the bedrock reaches a certain level,a wet bed is more likely to allow sliding and bed deformation (Irvine-Fynnet al.,2006).The formation of subglacial lakes at the terminus of some glaciers may be related to subglacier ablation caused by englacial meltwater (Nye,1960).

Figure 6 shows three large transverse crevasses on the east branch,from 4,600 m a.s.l.to 4,800 m a.s.l.(photos in Figure 9).Because the survey was done during the middle of the ablation period,there are many supraglacial water channels and considerable water import into the ice body through these surface crevasses.The form of temperate ice in this region is likely to be associated with infiltration of meltwater.

From the GPR data,we mapped the temperate ice thickness distribution (Figure 10).Temperate ice thickness at the mainstream line location was greater than that at the glacial edge,another indicator that meltwater gradually permeates into the bottom of the ice body under the influence of gravity.This suggests a wide range of hydrological connections.

Figure 7 Map of ice thickness (a) and bed landform (b) isoclines of L12 Glacier

Figure 8 Actual radar image (including an obvious temperate ice zone) (HZ segment in Figure 1,elevation between 4,390 and 4,580 m)

Figure 9 Photos of large traverse crevasses on the east branch.Distant view (a);Close shot (b)

Figure 10 Map of temperate ice thickness (unit: m)

5.4.Temperature and glacier change

To assess the impact of air temperature on L12 in recent years,we obtained annual mean air temperature (Figure 11)over the past 50 years from the Yemashan weather stations(39°30′N,96°30′E;4,250 m a.s.l.) near the glacier.Average temperatures from 1957 to 2008 have dramatically increased,with a big jump from 1997 to 1998.This trend continued after 2000.

Combined with previous investigations (Table 2),the rate of glacial area change from 1993 to 2000 has been twice as great as that from 1960 to 1990.This rate of change continued to increase after 2000.

It is generally believed that the movement of maritime glaciers is faster than continental glaciers,and continental glaciers move faster than extreme continental glaciers.These relative speeds have a close relationship with the water content of the ice body.We did not have enough temperature data from deep in the ice body to prove whether temperate ice is caused by temperature variation of ice bodies.As air temperatures rise,however,glacial surface meltwater will gradually increase,and the number of meltwater-importing ice bodies will gradually change englacial structures and hydrological characteristics.Meltwater will broaden englacial channels and further permeate the ice body.These processes will promote temperate ice formation.Our primary focus in the future will be on englacial temperatures,hydrological conditions,englacial pressure distributions and related themes.

Figure 11 Average air temperature over the past 50 years from Yemashan weather stations (39°30′N,96°30′E;4,250 m a.s.l.)

Table 2 Changes in glacial area for the L12 Glacier from 1960 to 2009

6.Conclusion

GPR is very effective at revealing complex subsurface structures in extreme continental glaciers.Englacial information such as crevasses,temperatures,ice thicknesses,and glacier beds was analyzed for the causes of temperate ice.It was concluded that:

(1) We used 2D FDTD to simulate certain typical englacial information.

(2) Through GPR surveys,we obtained the ice thickness and bed landform,which will be useful for prospective research.

(3) We determined that most temperate ice is distributed near the ice bed.Although the air temperature cannot directly affect the ice body,the infiltration of meltwater can change the thermal conditions of the ice body.

(4) The characteristics of the ice body are strongly affected by thickness,surface flow velocity,crevasses,meltwater,and air temperature.Each aspect all has some effects for glacial ablation.Therefore,in the future,we need combine each aspect to treat glacier ablation.

(5) Under the influence of crevasses,subglacial structure,and ice thickness,thickness reduction on the glacial east branch is more rapid than on the west branch.

Acknowledgments:

This work was funded by a project of the Chinese Glacier Inventory (41071047 and KZCX2-YW-GJ04).We express our thanks to Wu Wang,Gang Li,DongHui ShangGuan,and XiaoFeng Deng for guidance during the field survey.

REFERENCES

Alley R,Berntsen T,Bindoff NL,Chen ZL,Chidthaisong A,Friedlingstein P,Gregory J,Hegerl G,2007.Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change.In: Solomon SD,Qin M,Manning Z,Chen M,Marquis KB,Tignor AM,Miller HL (eds.).Climate Change 2007: The Physical Science Basis (Summary for Policymakers).Cambridge University Press,Cambridge,U.K.

Barry RG,2006.The status of research on glaciers and global glacier recession: A review.Progress in Physical Geography,30: 285–306.

Bartholomaus TC,Anderson RS,Anderson SP,2011.Growth and collapse of the distributed subglacial hydrologic system of Kennicott Glacier,Alaska,USA,and its effects on basal motion.Journal of Glaciology,57:985–1002.

Bergmann T,Robertsson JOA,Holliger K,1998.Finite-difference modeling of electromagnetic wave propagation in dispersive and attenuating media.Geophysics,63(3): 856–867.

Berthier E,2010.Volume loss from Bering Glacier (Alaska),1972–2003:Comment on Muskett and others (2009).Journal of Glaciology,56(197):555–557.

Binder D,Brückl E,Roch KH,Behm M,Schöner W,Hynek B,2009.Determination of total ice volume and ice-thickness distribution of two glaciers in the Hohe Tauern region,eastern Alps,from GPR data.Annals of Glaciology,50(51): 71–79.

Bingham RG,Nienow PW,Sharp MJ,Boon S,2005.Subglacial drainage processes at a high Arctic polythermal valley glacier.Journal of Glaciology,51(172): 15–24.

Björnsson H,1998.Hydrological characteristics of the drainage system beneath a surging glacier.Nature,395: 771–774.

Bongiovanni MV,Bonomo N,de la Vega M,Martino L,Osella A,2008.Rapid evaluation of multifrequency EMI data to characterize buried structures at a historical Jesuit mission in Argentina.Journal of Applied Geophysics,64(1–2): 37–46.

Bowersox R,2002.Hydrology of a Glacial Dominated System,Copper River,Alaska.University of California-Davis,pp.1–20.

Brandt O,Taurisano A,Giannopoulos A,Kohler J,2009.What can GPR tell us about cryoconite holes? 3D FDTD modeling,excavation and field GPR data.Cold Regions Science and Technology,55: 111–119.

Davis ME,Thompson LG,2006.An Andean ice-core record of a Middle Holocene mega-drought in North Africa and Asia.Annals of Glaciology,43: 34–41.

Du WT,Qin X,Liu YS,Wang XF,2008.Variation of the Laohugou Glacier No.12 in the Qilian Mountains.Journal of Glaciology and Geocryology,30(3): 373–379.

Fritzsche D,Schütt R,Meyer H,Miller H,Wilhelms F,Opel T,Savatyugin LM,2005.A 275-year ice-core record from Akademii Nauk ice cap,Severnaya Zemlya,Russian Arctic.Annals of Glaciology,42: 361–366.

Golledge NR,2007.Sedimentology,stratigraphy,and glacier dynamics,western Scottish highlands.Quaternary Research,68(1): 79–95.

Guan W,Hu HS,2008.Finite-difference modeling of the electroseismic logging in a fluid-saturated porous formation.Journal of Computational Physics,227(11): 5633–5648.

Haas C,Le Goff H,Audrain S,Perovich D,Haapala J,2011.Comparison of seasonal sea-ice thickness change in the Transpolar Drift observed by local ice mass-balance observations and floe-scale EM surveys.Annals of Glaciology,52(57): 97–102.

Hagg W,Braun L,2005.The influence of glacier retreat on water yield from high mountain areas: Comparison of Alps and Central Asia.In: de Jong C,Collins D,Ranzi R (eds.).Climate and Hydrology in Mountain Areas.John Wiley &Sons,Chichester,U.K.,pp.261–275.

Harper JT,Humphrey NF,Pfeffer WT,Fudge T,O’Neel S,2005.Evolution of subglacial water pressure along a glacier’s length.Annals of Glaciology,40(1): 31–36.

Herzfeld UC,Mayer H,2003.Seasonal comparison of ice-surface structures in the ablation area of Jakobshavn Isbrae drainage system,West Greenland.Annals of Glaciology,37: 199–206.

Hooke RLB,Fastook J,2007.Thermal conditions at the bed of the Laurentide ice sheet in Maine during deglaciation: Implications for esker formation.Journal of Glaciology,53(183): 646–658.

Humbert A,2010.The temperature regime of Fimbulisen,Antarctica.Annals of Glaciology,51(55): 56–64.

Iken A,Bindschadler RA,1986.Combined measurements of subglacial water pressure and surface velocity of Findelengletscher,Switzerland:Conclusions about drainage system and sliding mechanism.Journal of Glaciology,32(110): 101–119.

Irvine-Fynn TDL,Moorman BJ,Williams JLM,Walter FSA,2006.Seasonal changes in ground-penetrating radar signature observed at a polythermal glacier,Bylot Island,Canada.Earth Surface Processes and Landforms,31(7): 892–909.

Irvine-Fynn TDL,Moorman BJ,Willis IC,Sjogren DB,Hodson AJ,Mumford PN,Walter FSA,Williams JLM,2005.Geocryological processes linked to high Arctic proglacial stream suspended sediment dynamics:Examples from Bylot Island,Nunavut,and Spitsbergen,Svalbard.Hydrological Processes,19(1): 115–135.

Jørgensen AS,Andreasen F,2007.Mapping of permafrost surface using ground-penetrating radar at Kangerlussuaq Airport,western Greenland.Cold Regions Science and Technology,48(1): 64–72.

Joseph RM,Taflove A,1997.FDTD: Maxwell’s equations models for nonlinear electrodynamics and optics.Antennas and Propagation,45(3):364–374.

Khromova TE,Osipova GB,Tsvetkov DG,Dyurgerov MB,Barry RG,2006.Changes in glacier extent in the eastern Pamir,central Asia,determined from historical data and ASTER imagery.Remote Sensing of Environment,102(1–2): 24–32.

Kofman L,Ronen A,Frydman S,2006.Detection of model voids by identifying reverberation phenomena in GPR records.Journal of Applied Geophysics,59(4): 284–299.

Kotlyakov V,Macheret YY,1987.Radio echo-sounding of sub-polar glaciers in Svalbard: Some problems and results of Soviet studies.Annals of Glaciology,9: 151–159.

Lampe B,2003.Finite-difference time-domain modeling of ground-penetrating radar antenna systems.Ph.D.Thesis,Swiss Federal Institute of Technology,Zurich.

Liu CH,Kang ES,Liu SY,1999.Study on the glaciers variation and its runoff response in the arid region of northwest China.Science in China(Series D),29(supp.S1): 55–62.

Liu CH,Xie ZC,1988.Recent variation of glacier and forecast of its tendency in Qilian Mountains.Chinese Science Bulletin,33(8): 620–623.

Liu S,Xie ZC,Wang NL,Ye BS,1999.Mass balance sensitivity to climate change: A case study of glacier No.1 at Urumqi riverhead,Tianshan Mountains,China.Chinese Geographical Science,9(2): 134–140.

Liu SX,Zeng ZF,2007.FDTD simulation for ground penetrating radar wave in dispersive medium.Chinese Journal of Geophysics,50(1):299–306.

Macheret YY,Moskalevsky MY,Vasilenko EV,1993.Velocity of radio waves in glaciers as an indicator of their hydrothermal state,structure and regime.Journal of Glaciology,39(132): 373–384.

Molnia BF,2007.Late nineteenth to early twenty-first century behavior of Alaskan glaciers as indicators of changing regional climate.Global and Planetary Change,56(1–2): 23–56.

Moore PL,Iverson NR,Cohen D,2009.Ice flow across a warm-based/cold-based transition at a glacier margin.Annals of Glaciology,50(52): 1–8.

Murray T,Porter PR,2001.Basal conditions beneath a soft-bedded polythermal surge-type glacier: Bakaninbreen,Svalbard.Quaternary International,86(1): 103–116.

Nye JF,1960.The response of glaciers and ice-sheets to seasonal and climatic changes.Proceedings of the Royal Society (Series A): Mathematical,Physical &Engineering Sciences,256(1287): 559–584.

Pattyn F,Delcourt C,Samyn D,de Smedt B,Nolan M,2009.Bed properties and hydrological conditions underneath McCall Glacier,Alaska,USA.Annals of Glaciology,50(51): 80–84.

Paul F,Kääb A,Maisch M,Kellenberger T,Haeberli W,2004.Rapid disintegration of Alpine glaciers observed with satellite data.Geophysical Research Letters,31: L21402.

Pettersson R,2004.Dynamics of the cold surface layer of polythermal Storglaciären,Sweden.Ph.D.Thesis,Department of Physical Geography and Quaternary Geology,Stockholm University,Stockholm.

Pettersson R,Jansson P,Huwald H,Blatter H,2007.Spatial pattern and stability of the cold surface layer of Storglaciären,Sweden.Journal of Glaciology,53(180): 99–109.

Plewes LA,Hubbard B,2001.A review of the use of radio-echo sounding in glaciology.Progress in Physical Geography,25(2): 203–236.

Reid JR,Callender E,1965.Origin of debris-covered icebergs and mode of flow of ice into Miller Lake,Martin River Glacier,Alaska.Journal of Glaciology,5(40): 497–503.

Reineix A,Martinat B,Berthelier JJ,Ney R,2001.FDTD method for the theoretical analysis of the NETLANDER GPR.In: Conference on the Geophysical Detection of Subsurface Water on Mars.Houston,Texas,August 6–10,pp.88–89.

Rivera A,Casassa G,Bamber J,Kääb A,2005.Ice-elevation changes of Glaciar Chico,southern Patagonia,using ASTER DEMs,aerial photographs and GPS data.Journal of Glaciology,51(172): 105–112.

Rolstad C,Norland R,2009.Ground-based interferometric radar for velocity and calving-rate measurements of the tidewater glacier at Kronebreen,Svalbard.Annals of Glaciology,50(50): 47–54.

Serrano E,San José J,Agudo C,2006.Rock glacier dynamics in a marginal periglacial high mountain environment: Flow,movement (1991–2000)and structure of the Argualas Rock Glacier,the Pyrenees.Geomorphology,74(1–4): 285–296.

Shean DE,Marchant DR,2010.Seismic and GPR surveys of Mullins Glacier,McMurdo Dry Valleys,Antarctica: Ice thickness,internal structure and implications for surface ridge formation.Journal of Glaciology,56(195): 48–64.

Shi YF,Liu SY,2000.Estimation on the response of glaciers in China to the global warming in the 21st century.Chinese Science Bulletin,45(7):668–672.

Shugar DH,Rabus BT,Clague JJ,2010.Elevation changes (1949–1995) of Black Rapids Glacier,Alaska,derived from a multi-baseline InSAR DEM and historical maps.Journal of Glaciology,56(198): 625–634.

Su Z,Shi YF,2002.Response of monsoonal temperate glaciers to global warming since the Little Ice Age.Quaternary International,97–98:123–131.

Sun ZZ,Xie ZC,1981.The recent changes and trend of Glacier No.12 at the jokul of Qilian Mountains.Chinese Science Bulletin,26(6): 366–369.

Tamura T,Ohshima KI,Enomoto H,Tateyama K,Muto A,Ushio S,Massom RA,2006.Estimation of thin sea-ice thickness from NOAA AVHRR data in a polynya off the Wilkes Land coast,east Antarctica.Annals of Glaciology,44(1): 269–274.

Taurisano A,Tronstad S,Brandt O,Kohler J,2006.On the use of ground penetrating radar for detecting and reducing crevasse-hazard in Dronning Maud Land,Antarctica.Cold Regions Science and Technology,45(3): 166–177.

Thorsteinsson T,Waddington ED,Fletcher RC,2003.Spatial and temporal scales of anisotropic effects in ice-sheet flow.Annals of Glaciology,37(1): 40–48.

van der Veen CJ,Leftwich T,von Frese R,Csatho BM,Li J,2007.Subglacial topography and geothermal heat flux: Potential interactions with drainage of the Greenland ice sheet.Geophysical Research Letters,34:L12501.

Wang DM,McMechan GA,2002.Finite-difference modeling of borehole ground penetrating radar data.Journal of Applied Geophysics,49(3):111–127.

Wang ZT,Liu CH,You GX,Pu JC,Yang HA,Tian PY (eds.),1981.Glacier Inventory of China (I): Qilian Mountains.Lanzhou Institute of Glaciology and Cryopedology,Chinese Academy of Sciences,Lanzhou,China.

Welch BC,Jacobel RW,Arcone SA,2009.First results from radar profiles collected along the US-ITASE traverse from Taylor Dome to South Pole(2006–2008).Annals of Glaciology,50(51): 35–41.

Wohlleben T,Sharp M,Bush A,2009.Factors influencing the basal temperatures of a high Arctic polythermal glacier.Annals of Glaciology,50(52):9–16.

Wu Z,Liu SY,Zhang SQ,2009.Structural characteristics of the No.12 Glacier in Laohugou valley Qilian Mountain based on the ground penetrating radar sounding.Advances in Earth Science,24(10): 1149–1158.

October 17,2011 Accepted: March 26,2012

*Correspondence to: Zhen Wu,Key Laboratory of the Ecohydrology of Inland River Basins/Heihe Key Laboratory of Ecohydrology and Integrated Rivers.No.320,West Donggang Road,Lanzhou,Gansu 730000,China.Email: wuzhen@lzb.ac.cn