Contrasting thinning patterns between lake- and land-terminating glaciers in the Bhutan Himalaya

Despite the importance of glacial lake development in ice dynamics and glacier thinning, in situ and satellite-based measurements from lake-terminating glaciers are sparse in the Bhutan Himalaya, where a number of supraglacial lakes exist. We acquired in situ and satellite-based observations across lakeand land-terminating debris-covered glaciers in the Lunana 15 region, Bhutan Himalaya. A repeat differential global positioning system survey reveals that thinning of the debris-covered ablation area of the lake-terminating Lugge Glacier (−4.67 ± 0.02 m a) is more than three times greater than that of the landterminating Thorthormi Glacier (−1.40 ± 0.01 m a) for the 2004–2011 period. The surface flow velocities decrease downglacier along Thorthormi Glacier, whereas they increase from the upper part of the ablation area to the terminus of Lugge Glacier. Numerical experiments using a two-dimensional ice flow model demonstrate that the rapid thinning of Lugge Glacier 20 is driven primarily by a negative surface mass balance and that the dynamically induced change in ice thickness is small. However, the thinning of Thorthormi Glacier is suppressed by a longitudinally compressive flow regime. The magnitude of dynamic ice thickening compensates the glacier thinning, suggesting that over half of the negative surface mass balance is counterbalanced by the ice dynamics of Thorthormi Glacier. Multiple ponds on Thorthormi Glacier have been expanding since 2000 and merged into a single proglacial lake, with the glacier terminus detaching from its terminal moraine in 2011. Numerical 25 experiments suggest that the speed up and thinning of Thorthormi Glacier will accelerate with continued proglacial lake development.


Introduction 28
The spatially heterogeneous shrinkage of Himalayan glaciers has been revealed by in situ measurements (Yao et  2013). Bhutanese glaciers are inferred to be particularly sensitive to changes in air temperature and precipitation 37 because they are affected by monsoon-influenced, humid climate conditions (Fujita and Ageta, 2000;Fujita, 2008; 38 Sakai and Fujita, 2017). Mass loss of Gangju La Glacier in central Bhutan was much greater than those of glaciers 39 in the eastern Himalaya and southeastern Tibet for the recent decade (Tshering and Fujita, 2016). It is crucial to 40 investigate the mechanisms driving the mass loss of Bhutanese glaciers to provide more insight for the glacier mass 41 balance (Zemp et al., 2015) and improve projections of global sea level rise and glacier evolution (Huss and Hock, 42 2018). 43 In recent decades, glacial lakes have formed and expanded at the termini of retreating glaciers in the Himalayas 44 Widespread thinning of Himalayan glaciers has been revealed by differencing multitemporal DEMs constructed 69 from satellite image photogrammetry because the surface of debris-covered glaciers can be highly variable, making 70 access difficult to get large amounts of data (e.g., Gardelle et al., 2013;Maurer et al., 2016;Brun et al., 2017). In 71 particular, unmanned autonomous vehicle (UAVs) is a powerful tool to obtain higher-resolution imagery than 72 satellite, and thus resolves the highly variable topography and thinning rates of debris-covered surface more 73 accurately (e.g., Immerzeel et al., 2014;Vincent et al., 2016). Repeated differential GPS (DGPS) measurements, 74 which are acquired with centimetre-scale accuracy, also enable us to evaluate elevation changes of several metres 75 (e.g., Fujita et al., 2008). Although their temporal and spatial coverage can be limited, the repeated DGPS 76 measurements have been successfully acquired to investigate the surface elevation changes of debris-free glaciers 77 in Bhutan (Tshering and Fujita, 2016) and the Inner Tien Shan . 78 This study aims to quantify the contributions of ice dynamics and SMB to the thinning of adjacent land-and 79 lake-terminating glaciers. To investigate the importance of glacial lake formation and expansion on glacier thinning, 80 we measured surface elevation changes on a lake-and a land-terminating glacier in the Lunana region, Bhutan 81 9 thermal resistance (e.g., Reid and Brock, 2010). Using eight ASTER images (90 m resolution, Level 3A1 data) 215 obtained between October 2002 and October 2010 (Table S4), along with the NCEP/NCAR reanalysis climate data 216 (NCEP-2, Kanamitsu et al., 2002), we calculated the distribution of mean thermal resistance on the two target 217 glaciers. Surface albedo is calculated using three-visible near-infrared sensors (VNIR; bands 1-3), and surface 218 temperature is obtained from an average of five sensors in the thermal infrared (TIR; bands 10-14). Automatic 219 weather station (AWS) observations from the terminal moraine of Lugge Glacial Lake (4524 m a.s.l., Fig. 1a) 220 showed that the annual mean air temperature during 2002-2004 was ∼0 °C, and annual precipitation was 900 mm 221 in 2003 (Suzuki et al., 2007b). The air temperature at the AWS elevation was estimated using the pressure level 222 atmospheric temperature and geopotential height (Sakai et al., 2015), and then modified for each 90 × 90 m mesh 223 grid points using a single temperature lapse rate (0.006 ºC km −1 ). The wind speed was assumed to be 2.0 m d −1 , 224 which is the two-years average of the 2002-2004 AWS record (Suzuki et al., 2007b). The uncertainties in the 225 thermal resistance and albedo were evaluated as 107 and 40%, respectively, by taking the standard deviations 226 calculated from multiple images at the same location (Fig. S2). 227 The SMB of the debris-covered ablation area was calculated by a heat and mass balance model that included 228 debris-covered effects (Fujita and Sakai, 2014). First, the surface temperature is determined to satisfy Eq. (2) using 229 the estimated thermal resistance and an iterative calculation, and then, if the heat flux toward the ice-debris 230 interface is positive, the daily amount of ice melt beneath the debris mantle ( 3 , kg m −2 d −1 ) is obtained as follows: 231 where @ is the length of a day in seconds (86400 s), and B is the latent heat of fusion of ice (3.33 × 105 J kg -1 ). 234 Annual mass balance of debris-covered part ( , m w.e. a -1 ) is expressed as: here 5 and G are snow and rain, respectively, which are distinguished from precipitation depending on air 238 temperature. Evaporation from debris and snow surfaces is expressed in the same formula but they are calculated in 239 10 different schemes because temperature and saturation conditions of the debris and snow surfaces are different. 3 240 and 5 are the daily discharge from the debris and snow surfaces, respectively. Discharge and evaporation from the 241 snow surface was calculated only when snow layer was formed on the debris. Because snow layer does not exist at 242 the end of melting season in the current climate condition and at the elevation of debris-covered area, snow 243 accumulation ( 5 ) is compensated with evaporation and discharge from snow surface during a calculation year. show no clear anomaly against the long-term mean SMB (1979-2017) (Fig. S4). 260 The sensitivity of the simulated meltwater was evaluated against the meteorological parameters used in the SMB 261 model. We chose meltwater instead of SMB to quantify the uncertainty because the SMB uncertainty cannot be 262 evaluated as absolute value. The tested parameters are surface albedo, air temperature, precipitation, relative 263 humidity, solar radiation, thermal resistance and wind speed. The thermal resistance and albedo uncertainties were 264 based on the standard deviations derived from the eight ASTER images used to estimate these parameters (Fig. S2). 265 Each meteorological variable uncertainty, with the exceptions of the thermal resistance and albedo uncertainties, 266 was assumed to be the root mean square error (RMSE) of the ERA-Interim reanalysis data against the observational 267 data (Fig. S3). The simulated meltwater uncertainty was estimated as the variation in meltwater within a possible 268 parameter range via a quadratic sum of the results from each meteorological parameter. 269 where ̇1^ and 1^ are the components of the strain rate and deviatoric stress tensors, respectively, and # is the 289 effective stress, which is described as 290

Ice dynamics
292 The rate factor ( ) and flow law exponent ( ) are material parameters. We used the commonly accepted value of 293 = 3 for the flow law exponent and employed a rate factor of = 75 MPa fO a fS , which was previously used to 294 model a temperate valley glacier (Gudmundsson, 1999). We assumed the glaciers were temperate because there 295 was no available information on the thermal states of the studied glaciers. The upper surface of the domain was assumed to be stress free, through which the ice flux was prescribed from 313 the surface velocity obtained by the satellite analysis. We assume no basal sliding and quadratic function (4th 314 order) for the velocity profile from the surface to the bed. The basal sliding velocity ( J ) was given as a linear 315 function of the basal shear traction ( XY,J ): 316 318 where is the sliding coefficient. We used constant sliding coefficients of C = 766 and 125 m a −1 MPa −1 over the 319 entire domains of Thorthormi and Lugge glaciers, respectively. These parameters were obtained by minimising the 320 RMSE between the modelled and measured surface flow velocities over the entire model domains (Fig. S5). 321

Experimental configurations 322
To quantify the effect of glacier dynamics on ice thickness change, we performed two experiments for Thorthormi 323 and Lugge glaciers. Experiment 1 was performed to compute the ice flow velocity fields under the present terminus 324 conditions. In this experiment, Thorthormi Glacier was treated as a land-terminating glacier by prescribing zero 325 horizontal velocity at the glacier front, whereas Lugge Glacier was treated as a lake-terminating glacier by applying 326 hydrostatic pressure at the front as a function of water depth. A stress-free boundary condition was given to the 327 calving front above the lake level. We employed glacier surface elevation in 2001 and water level of supraglacial 328 ponds and proglacial lake observed in 2004 as boundary conditions . 329 Experiment 2 was designed to investigate the influence of proglacial lakes on glacier dynamics. For Thorthormi 330 Glacier, we simulated a calving front with thickness of 106 m. Position of the hypothetical calving front was 331 determined at the place where only one lake depth was acquired from bathymetry survey in September 2011. The 332 surface level of the proglacial lake was assumed to be 4432 m a.s.l., which is the mean surface level of the 333 supraglacial ponds measured in September 2004 . Hydrostatic pressure and stress-free 334 conditions were applied to the lower boundary below and above the lake level, respectively. For Lugge Glacier, we 335 simulated a lake-free situation, with ice flowing to the contemporary terminal moraine, so that the glacier 336 terminates on land. Bedrock topography is derived from the bathymetric map (white lines in Fig. 1b Glacier, the flow velocities decrease down-glacier, ranging from ~110 m a −1 at the foot of the icefall to < 10 m a −1 381 at the terminus (Fig. 3a). The flow velocities of Lugge Glacier increase down-glacier, ranging from 20-60 to 50-80 382 m a −1 within 2000 m of the calving front (Fig. 3b). The flow velocity uncertainty was estimated to be 12.1 m a −1 , as 383 given by the mean off-glacier displacement from 3 February 2006 to 30 January 2007 (362 days) (Fig. S9). 384

Changes in glacial lake area 385
The supraglacial pond area near the front of Thorthormi Glacier progressively increased from 2000 to 2017, at a 386 mean rate of 0.09 km 2 a −1 while Lugge Glacial Lake also expanded from 2000 to 2017, at a mean rate of 0.03 km 2 387 a −1 (Fig. 4). The total area changes from 2000 to 2017 are 1.79 km 2 and 0.46 km 2 for Thorthormi and Lugge 388 glaciers, respectively. 389

Mass balance of the debris-covered surface 390
The simulated SMBs over the ablation area were −7.36 ± 0.12 m w.e. a −1 for Thorthormi Glacier and −5.25 ± 0.13 391 m w.e. a −1 for Lugge Glacier (Fig. 1c, Table 1). The debris-free surface has a more negative SMB than the debris-392 covered regions of the glaciers. The mean SMBs of the debris-free and debris-covered surfaces in the ablation area 393 of Thorthormi Glacier are −9.31 ± 0.68 and −7.30 ± 0.13 m w.e. a −1 , respectively, while those of Lugge Glacier are 394 −7.33 ± 0.41 m w.e. a −1 and −5.41 ± 0.18 m w.e. a −1 , respectively (Table 1). The sensitivity of simulated meltwater 395 in the SMB model was evaluated as a function of the RMSE of each meteorological variable across the debris-396 covered area (Fig. S10). Ice melting is more sensitive to solar radiation and thermal resistance. The influence of 397 thermal resistance on meltwater formation is considered to be small since the debris cover is sparse over the 398 glaciers. The estimated meltwater uncertainty is < 50% across most of Thorthormi and Lugge glaciers (Fig. S11). 399

Numerical experiments of ice dynamics 400
The ice thinning of Lugge Glacier was three times faster than that of Thorthormi Glacier. However, the mean SMB 401 was 1.4 times more negative at Thorthormi Glacier, suggesting a substantial influence of glacier dynamics on ice 402 thickness change. To quantify the contribution of ice dynamics to the ice thickness change, we performed 403 numerical experiments with the present (Experiment 1) and prescribed (Experiment 2) glacier geometries. 404

Experiment 1 -present terminus conditions 405
Modelled results for the present geometry show significantly different flow velocity fields for Thorthormi and 406 Lugge glaciers (Figs. 5c and 5d). Thorthormi Glacier flows faster (> 150 m a −1 ) in the upper reaches, where the 407 surface is steeper than the other regions (Fig. 5c). Down-glacier of the icefall, where the glacier surface is flatter, 408 the ice motion slows in the down-glacier direction, with the flow velocities decreasing to < 10 m a −1 near the 409 terminus (Fig. 5e). Ice flows upward relative to the surface across most of the modelled region (Fig. 5c). In contrast 410 to the down-glacier decrease in the flow velocities at Thorthormi, the computed velocities of Lugge Glacier are up 411 to ~40 m a −1 within 500-2200 m of the terminus, and it sharply increases to ~80 m a −1 at the calving front (Fig. 5f). 412 Ice flow is nearly parallel to the glacier surface, except for the more downward motion near the calving front (Fig.  413   5d). Within 3000 m of the terminus of Thorthormi Glacier, the modelled surface flow velocities are in good 414 agreement with the satellite-derived flow velocities (Fig. 5e). The calculated surface flow velocities of Lugge 415 Glacier agreed with the satellite-derived flow velocities within 16 % (Fig. 5f). with Experiment 1 (Figs. 5d and 6d). The upward ice motion appears within 3000 m of the terminus. The numerical 424 experiments demonstrate that the formation of a proglacial lake causes significant changes in ice dynamics. 425

Simulated surface flow velocity uncertainty 426
Basal sliding accounts for 97 % and 75 % of the simulated ice flow velocity in the ablation area of Thorthormi 427 and Lugge glaciers, respectively (Figs. 5e and 5f), suggesting that ice deformation is insufficient to represent the ice 428 flow regardless of assumption of ice temperature. Standard deviations of ASTER-derived surface velocities are 2.9 429 and 6.7 m a −1 for Thorthormi and Lugge glaciers, respectively, which are considered as interannual variability in 430 the measured surface velocities (Fig. 3). The RMSEs between the modelled and measured flow velocities were 431 computed as a measure of the model performance. For Thorthormi Glacier, the model exhibits similar sensitivities 432 to the sliding coefficient and ice thickness while the model is more sensitive to the ice thickness than the sliding 433 coefficient for Lugge Glacier (Fig. S5). Uncertainty of the emergence velocity is affected by factors such as 434 accuracy of flow velocity measurement, interannual variability, and RMSE between modelled and measured flow 435 velocity so that we performed sensitivity test by changing ±30 % of ice thickness and sliding coefficient. It shows 436 that the simulated surface flow velocities of Thorthormi Glacier vary by ±30 % when the constant sliding 437 coefficient (C) and ice thickness are varied by ±30 % (Fig. S12) (Fig. 7a). 445 In contrast, the emergence velocity of Lugge Glacier was −1.69 ± 0.18 m a −1 within 700-2500 m of the terminus 446 ( Fig. 7a) and greatly negative value at the calving front (−10 m a −1 ) due to the increase in surface flow velocities 447 toward the glacier front (Fig. 5f). Under the Experiment 1 conditions, the estimated ℎ⁄ are −4.88 ± 0.34 m a −1 448 within 4200 m of Thorthormi Glacier and −7.46 ± 0.32 m a −1 within 700-2500 m from the calving front of Lugge 449 Glacier (Fig. 7). 450 The emergence velocity computed under contrasting geometries (Experiment 2) varies from that with the present 451 geometries (Experiment 1) for both Thorthormi and Lugge glaciers. For the lake-terminating condition of 452 Thorthormi Glacier, the mean emergence velocity turns into negative (−1.37 ± 0.23 m a −1 ) within 3700 m of the 453 terminus. The mean emergence velocity of Lugge Glacier computed with the land-terminating condition is less 454 negative (−0.78 ± 0.28 m a −1 ) within 700-2500 m of the terminus. Given the same SMB distribution, the mean 455 ℎ⁄ is computed as −9.46 ± 0.36 m a −1 for Thorthormi Glacier with the lake-terminating condition and −6.55 ± 456 0.42 m a −1 for the land-terminating Lugge Glacier (Table 1). 457 parts of nine lake-terminating glaciers in the Everest area (approximately −2.5 m a −1 ), which was 67% more 476 negative than that of 18 land-terminating glaciers (approximately −1.5 m a −1 ). The ℎ⁄ of lake-terminating 477 glaciers in Yakutat ice field, Alaska (−4.76 m a −1 ) was ~30% more negative than the neighbouring land-terminating 478 glaciers (Trüssel et al., 2013). It should be noted that the difference in ℎ⁄ between Lugge and Thorthormi 479 glaciers derived from DGPS-DEMs (3.3 times) is much greater than the numbers previously reported in the Nepal 480 Himalaya, suggesting that ice dynamics play a more significant role here. 481 20

Influence of ice dynamics on glacier thinning 482
The estimated ℎ⁄ are 3.5 time more negative than the DPGS observation for Thorthormi Glacier and 60 % 483 more negative than the DGPS observation for Lugge Glacier, respectively (Table 1) Glacier is 40 % more negative than that of Lugge Glacier. Since debris cover looks sparse across the ablation area 489 of both glaciers (Fig. S1), the more negative SMB of Thorthormi Glacier could be explained by the glacier situated 490 at lower elevations (Fig. 2b). The calculated SMBs (Thorthomi < Lugge) and observed ℎ⁄ (Lugge < 491 Thorthormi) suggest that contribution of glacier dynamics is substantially different in the two glaciers. The 492 horizontal flow velocities of Lugge Glacier are nearly uniform along the central flowline, with ice flow parallel to 493 the glacier surface (Fig. 5d), suggesting that the dynamically induced ice thickness change is small. The computed 494 emergence velocity is negative (−1.69 ± 0.18 m a −1 ), which means the ice dynamics rather accelerates glacier 495 thinning. In contrast to Lugge, the flow velocities of Thorthormi Glacier decrease toward the terminus (Fig. 5c) Thorthormi Glacier is equivalent to 60 % of the negative SMB, implying that one-third of the surface ablation is 499 counterbalanced by ice dynamics. In other words, dynamically induced ice thickening partly compensates for the 500 negative SMB. 501 Experiment 1 demonstrates that the difference in emergence velocity between the land-and lake-terminating 502 glaciers leads to contrasting thinning patterns. On the other hand, Experiment 2 demonstrates that the emergence 503 velocity was less negative (−0.78 ± 0.28 m a −1 ) in the absence of a glacial lake in Lugge Glacier, resulting in a 504 decrease in the thinning rate by 12 % as compared with the lake-terminating condition. The negative emergence 505 velocity suggests that Lugge Glacial Lake could have been inevitably formed, and the more negative emergence 506 velocity caused by the development of the lake would have accelerated the thinning of Lugge Glacier. For 507 Thorthormi Glacier, the emergence velocity under the lake-terminating condition is negative (−1.37 ± 0.23 m a −1 ), 508 resulting in a doubled thinning rate (4.88 to 9.46 m a −1 , Table 1). Our ice flow modelling demonstrates that thinning 509 21 will be accelerated in association with the development of a supraglacial lake in the terminal part of Thorthormi 510 Glacier. 511 Contrasting patterns of glacier thinning and horizontal flow velocities between the land-and lake-terminating 512 glaciers are consistent with satellite-based observations over lake or ocean-terminating glaciers and neighbouring 513 land-terminating glaciers in the Nepal Himalaya (King et al., 2017) and Greenland (Tsutaki et al., 2016). A 514 decrease in the down-glacier flow velocities over the lower reaches of land-terminating glaciers suggests a 515 longitudinally compressive flow regime, which would result in a positive emergence velocity and therefore 516 thickening to compensate for the negative SMB. Conversely, for lake-terminating glaciers, an increase in the down-517 glacier flow velocities suggests a longitudinally stretching flow regime, which would yield a negative emergence 518 velocity, resulting in accelerated ice thinning. Contrasting flow regimes modelled in this study suggest that the 519 mechanisms would not only be applicable to Thorthormi and Lugge glaciers, but to other lake-and land-520 terminating glaciers worldwide where contrasting thinning patterns are observed. 521 The thinning rate calculated from the model is ~5 m a −1 more negative than the observation over the entire 522 domain of Lugge Glacier and also the lower part of Thorthormi Glacier (Fig. 7b), which is probably due to the 523 uncertainties in the estimated ice thickness and basal sliding conditions and SMB. The two-dimensional feature is 524 another reason for the insufficient modelled results because the model neglects drag from the side walls and 525 changes in glacier width. The SMB uncertainty is < 50% over a large portion of Thorthormi and Lugge glaciers 526 (Fig. S11). Nevertheless, our numerical experiments suggest that dynamically induced ice thickening compensates 527 the negative SMB in the lower part of land-terminating glacier, resulting in less ice thinning in contrast to the lake-528 terminating glacier. Further measurements of the spatial distributions of ice thickness and SMB will help in 529 deriving more accurate estimates of the effect of ice dynamics on glacier thinning. 530

Proglacial lake development and glacier retreat 531
Lugge Glacial Lake has expanded continuously and at a nearly constant rate from 2000 to 2017 (Fig. 4). 532 Bathymetric data suggest that glacier ice below the lake level accounted for 89% of the full ice thickness at the 533 calving front in 2002 (Fig. 5b). If lake level is close to the ice flotation level, where the basal water pressure equals 534 the ice overburden pressure, calving caused by ice flotation regulates glacier front position (van der Veen, 1996), 535 and glacier could rapidly retreat (e.g., Motyka (Fig. 3a), it is likely that the thinning and retreat of Thorthormi Glacier will be 561 accelerated in the near future due to the formation and expansion of the proglacial lake. 562

Conclusions 563
To better understand the importance of glacial lake formation on rapid glacier thinning, we carried out field and 564 satellite-based measurements across the lake-terminating Lugge Glacier and the land-terminating Thorthormi 565 Glacier in the Lunana region, Bhutan Himalaya. Surface elevations were surveyed in 2011 by differential GPS 566 (DGPS) across the lower parts of the glaciers and compared with a 2004 DGPS survey. Changes in surface 567 elevation were also measured by differencing satellite-based DEMs. The flow velocity and area of glacial lake were 568 determined from optical satellite images. We also performed numerical experiments to quantify the contributions of 569 surface mass balance (SMB) and ice dynamics in relation to the observed ice thinning. 570 Lugge Glacier has experienced rapid ice thinning which is 3.3 times greater than Thorthormi Glacier, even 571 though the SMB was less negative. The numerical modelling results, using the present glacier geometries, 572 demonstrate that Thorthormi Glacier is subjected to a longitudinally compressive flow regime, suggesting that 573 dynamically induced vertical extension compensates for the negative SMB, and thus results in less ice thinning than 574 at Lugge Glacier. Conversely, the computed negative emergence velocity suggests that the rapid thinning of Lugge 575 Glacier was driven by both surface melt and ice dynamics. This study reveals that contrasting ice flow regimes 576 cause different ice thinning observations between the lake-and land-terminating glaciers in the Bhutan Himalaya.