1. Introduction
The High Mountain Asia (HMA) region, known as the ‘Third Pole’, contains the largest concentration of ice outside of the polar ice sheets. The glaciers of HMA feed major water systems, which provide water and sanitation for over a billion people (Scott and others, Reference Scott, Zhang, Mukherji, Immerzeel, Mustafa, Bharati, Wester, Mishra, Mukherji and Shrestha2019). In particular, the Karakoram is the most heavily glaciated mountain range in Asia (RGI Consortium, 2017) and is a critical water source for large parts of Pakistan and parts of northern India (Scott and others, Reference Scott, Zhang, Mukherji, Immerzeel, Mustafa, Bharati, Wester, Mishra, Mukherji and Shrestha2019). However, climate change has led to increasingly negative mass balance, putting the area’s future at risk (Rounce and others, Reference Rounce, Hock and Shean2020; Zhang and others, Reference Zhang2023). Glacial lake outburst floods (GLOFs) in the region have also caused significant loss of human lives and infrastructure damage in recent decades (Shrestha and others, Reference Shrestha2023), and the risk of exposure to local communities and infrastructure due to growing proglacial lakes may potentially increase (Zheng and others, Reference Zheng2021). GLOFs in the Karakoram region occur through breaches of moraine or ice dams, which are associated with rapid (re)-organization of subglacial waters and channels (Nye, Reference Nye1976; Gudmundsson and others, Reference Gudmundsson, Björnsson and Pálsson1995; Flowers and others, Reference Flowers, Björnsson, Pálsson and Clarke2004; Kingslake and Ng, Reference Kingslake and Ng2013). Proglacial and proximal (ice-dammed) lakes, which are often hydraulically connected with the subglacial drainage network, also exert an important boundary condition for subglacial hydrology (Armstrong and Anderson, Reference Armstrong and Anderson2020).
The Karakoram region is also home to a high concentration of surge-type glaciers (Copland and others, Reference Copland2011; Sevestre and Benn, Reference Sevestre and Benn2015). Surges are a phenomenon characterized by cyclical, order-of-magnitude accelerations of glaciers that can be sustained for months to years (Meier and Post, Reference Meier and Post1969). They occur in geographical clusters that fall in ‘climatic envelopes’ that may provide favorable conditions for surge motion (Sevestre and Benn, Reference Sevestre and Benn2015). Buildups of basal water pressure are thought to play a role in the initiation and sustenance of surge motion (Kamb, Reference Kamb1987; Björnsson, Reference Björnsson1998; Flowers and others, Reference Flowers, Roux, Pimentel and Schoof2011). Surging is also associated with till deformation—e.g. Truffer and others (Reference Truffer, Harrison and Echelmeyer2000) and Minchew and Meyer (Reference Minchew and Meyer2020). However, the causes of surge behavior remain unclear, as not all surging glaciers seem to be directly attributable to changes in mass-balance state or thermal regime (Murray and others, Reference Murray2000; Liu and others, Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024).
Subglacial hydrology controls ice velocity through changes in effective pressure, defined as the difference between the ice overburden pressure and the water pressure at the bed (Cuffey and Paterson, Reference Cuffey and Paterson2010). Seasonal variations in subglacial hydrology modulate ice sheet and glacier velocities (Iken and others, Reference Iken, Röthlisberger, Flotron and Haeberli1983; Schoof, Reference Schoof2010). Numerous studies have shown that the velocity of glaciers is influenced by seasonal cycles (Zwally and others, Reference Zwally, Abdalati, Herring, Larsen, Saba and Steffen2002; Hart and others, Reference Hart, Young, Baurley, Robson and Martinez2022; Nanni and others, Reference Nanni, Scherler, Ayoub, Millan, Herman and Avouac2023). In alpine glaciers of HMA, observed regional speedups have been proposed to occur due to changes in melt volume and subglacial drainage efficiency. In particular, these glaciers can also exhibit a pattern of speedups in both the spring and fall (Beaud and others, Reference Beaud, Aati, Delaney, Adhikari and Avouac2022, Nanni and others, Reference Nanni, Scherler, Ayoub, Millan, Herman and Avouac2023). These studies suggest that seasonal speedups occur due to increases in meltwater production and subsequent lubrication at the ice–bed interface.
While surges and outburst flooding have, for the most part, been investigated as separate phenomena, multiple studies in the Karakoram have observed GLOFs to occur concurrently with transitions in surge motion, suggesting that subglacial hydrology may play a role in the synchronous timing of these events (Round and others, Reference Round, Leinss, Huss, Haemmig and Hajnsek2017; Bazai and others, Reference Bazai2022; Beaud and others, Reference Beaud, Aati, Delaney, Adhikari and Avouac2022). Understanding the role that subglacial hydrology plays in the severity and timing of these hazards could improve early warning systems for water availability and outburst flooding. While several in situ observational studies have been conducted and are in progress (e.g. Gilbert and others, Reference Gilbert2020; Pritchard and others, Reference Pritchard, King, Goodger, McCarthy, Mayer and Kayastha2020; Miles and others, Reference Miles, Hubbard, Quincey, Miles, Irvine-Fynn and Rowan2019), there are very few direct observations of subglacial hydrology in HMA. Therefore, in this study, we lay the groundwork for investigating the role of subglacial hydrology in ice dynamics and outburst flooding through modeling.
We focus on Shishper Glacier (36.40∘N 74.61∘E) in the eastern Karakoram range in Pakistan (Fig. 1). The glacier has also been referred to in literature as Shisper and Shishpare. Located in the Hunza Valley in Gilgit-Baltistan, Pakistan, Shishper is part of a surge and lake drainage system with another glacier to its west, called Muchuwar (also previously spelled as Muchuhar or Mochowar). The two glaciers were connected prior to 1950, when the two separated (Muhammad and others, Reference Muhammad2021). Shishper’s main trunk is 7 km long and is fed by several tributary glaciers at the northeast (upper-elevation) side. In total, the glacier is 15 km in length.
Figure 1. (a) Outlines of adjacent valley glaciers Shishper and Muchuwar (Randolph Glacier Inventory Version 6.0) overlaid on Landsat 8 OLI NIR imagery from December 2016. Our modeled domain is outlined in red. (b) Surface elevation from TanDEM-X 90m DEM, with contours showing the terrain elevation in meters. (c) Ice thickness from the global dataset by Millan and others Reference Millan, Mouginot, Rabatel and Morlighem(2022). All coordinates are projected to WGS 84/UTM Zone 42N.
Both Shishper and Muchuwar have surged cyclically for as long as observations have been recorded, since the early 1900s (Beaud and others, Reference Beaud, Aati, Delaney, Adhikari and Avouac2022). Shishper underwent major surges in 1973, 2000–2011 and most recently between 2017-19 (Bhambri and others, Reference Bhambri2020). During this time, the terminus advanced 1.5 km (Bhambri and others, Reference Bhambri2020). In June 2019, the surge and subsequent lake drainage resulted in the closing of two power plants, the evacuation and considerable damage to some houses in the downstream village and lasting damage to agricultural land, and finally, the destruction of the main road bridge crossing the stream, affecting transport along the main transport axis in the region. In mid-November 2018, the advancement of Shishper blocked meltwater flow from Muchuwar Glacier, which created an ice-dammed proximal lake (Beaud and others, Reference Beaud, Aati, Delaney, Adhikari and Avouac2022). This lake tends to fill up in November–December and in May to a depth of 30–80m, with an estimated volume of 30 million m3. When the lake drains, the outburst flood drains through Shishper’s terminus and down into the valley below. The maximum river flow observed at the downstream village of Hassanabad is 150–200 m3 s−1, compared to a base flow of 20 m3 s−1 (Muhammad and others, Reference Muhammad2021). After the lake is filled in the winter, drainage occurs more gradually, as opposed to the spring filling, which results in a more dramatic drainage of the lake.
In this study, we simulate the seasonal dynamics of the subglacial drainage system of Shishper Glacier in isolation from velocity coupling and lake drainage. We use a state-of-the-art subglacial hydrology model, forced with realistic meltwater inputs, to gain insight into the evolution of the water flow and pressure distribution beneath the glacier. The following sections describe the modeling methods and assumptions, meltwater forcing data, simulation results, and limitations of the approach.
2. Model setup and assumptions
To simulate the subglacial hydrological system of Shishper Glacier, we employ the SHAKTI (Subglacial Hydrology and Kinetic, Transient Interactions) model (Sommers and others, Reference Sommers, Rajaram and Morlighem2018), which is implemented in the Ice-sheet and Sea-level System Model (ISSM) (Larour and others, Reference Larour, Seroussi, Morlighem and Rignot2012). The current implementation of SHAKTI in ISSM is the simplified formulation described in Sommers and others Reference Sommers2023, which neglects englacial storage, opening by sliding, and melt due to changes in the pressure melting point. SHAKTI is capable of modeling a variety of network systems between the end-member cases of efficient and inefficient drainage systems. It does this by allowing the hydraulic transmissivity to vary spatially and temporally. In addition, it accounts for varying laminar, turbulent, and intermediate flow regimes.
The model domain is traced from the Randolph Glacier Inventory, Version 6.0 (RGI Consortium, 2017). The tributary branches of Shishper Glacier, located above 3500 m a.s.l., likely experience less liquid precipitation and decreased melting compared to the lower section of the main trunk and therefore may not contribute significantly to the subglacial hydrological system. Our aim is to examine the evolution in the hydrology in the main trunk, rather than evaluating the exact quantity of subglacial water in the system; for these reasons, we reserve including hydrological contributions from the tributary glaciers for future work. Furthermore, we neglect frictional heating due to basal sliding, which may decrease the flux of meltwater through the hydrological system; however, our intention is to isolate the effects of seasonal melt on the drainage system, so we also reserve calculation of melt from frictional sliding for future work. The modeled hydrological domain overlaid on the RGI 6.0 outline is shown in Fig. 1, depicting a glacier outline from 2016, before the 2017-19 surge event. We focus on modeling the subglacial hydrology for a steady geometry, with the glacier outline and ice thickness held constant throughout the simulations.
To obtain surface and bed geometries for the glacier, we use the TanDEM-X global DEM (German Aerospace Center, 2018) along with a global glacier thickness dataset (Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022). Glacier thickness is subtracted from surface elevation to obtain a bed topography, and all spatial data are projected to WGS 84/UTM Zone 42N. Radar mapping of subglacial topographies and glacier thicknesses at other glaciers have revealed large uncertainties associated with this ice thickness dataset, which was calculated using mass conservation techniques (Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022; Tober and others, Reference Tober, Christoffersen, Holt, Truffer and Larsen2024). In addition, it is likely that artificially smooth bed topographies calculated from mass conservation inversions may affect the results of subglacial hydrology simulations (MacKie and others, Reference MacKie, Schroeder, Zuo, Yin and Caers2021). Due to the lack of in situ observations to validate the model by Millan and others Reference Millan, Mouginot, Rabatel and Morlighem(2022), we emphasize that the exact routing of subglacial channels in our simulations is subject to the uncertainty associated with the estimated bed topography.
We manually trace the model domain to the RGI outline using the in-built functionality in ISSM. The DEM and bed topography data are interpolated onto a two-dimensional unstructured triangular mesh with 40 m resolution. The ideal mesh size and geometry were determined after conducting a winter equilibration for 600 days at varying mesh sizes (shown in Appendix B). We conclude from these tests that the location of channel formations is insensitive to mesh size. The 40 m resolution, which yields a mesh containing 3302 vertices and 6035 elements, provides enough detail and numerical stability while saving on computational costs. The mesh provides the basis for the P1 triangular Lagrange finite element solver used by SHAKTI. We test 20 slightly varying domain outlines with slightly different variations in domain outline, conducting a winter equilibration wherein all subglacial water is generated by basal melt (see Section ‘Establishing winter base state’) for 1000 days on each. The final geometry used for the transient simulations is chosen based on the criteria that mean ice–bed gap height, gap-integrated basal water flux (the approximate momentum equation for water velocity integrated over the gap height—see Sommers and others, Reference Sommers, Rajaram and Morlighem2018), and effective pressure equilibrate after 1000 days without anomalous numerical artefacts near corners or curvatures. The ice velocity is set to 0 throughout all of the transient simulations, isolating the seasonal evolution of subglacial hydrology without frictional heating feedbacks from basal sliding. All simulations in this work are carried out with ISSM Version 4.23 using a MATLAB interface on MacOS.
2.1. Surface melt timeseries
To estimate the timing and magnitude of seasonal meltwater inputs to the bed, we use the European Centre for Medium-Range Weather Forecasts (ECMWF)’s Reanalysis v5 (ERA5) (Muñoz-Sabater and others, Reference Muñoz-Sabater2021) as inputs to the temperature-indexed model by Litt and others Reference Litt(2019) to obtain spatio-temporally varying estimates for surface melt across the domain (Fig. 2). These ERA5 weather data are based on an array of field stations and weather models (Hersbach and others, Reference Hersbach2020) and directly provide estimates for snow cover, air temperature and total liquid precipitation across the 5 years (Muñoz-Sabater and others, Reference Muñoz-Sabater2021). Ice melt across the mesh is calculated using the temperature index (TI) melt parametrization from Litt and others (Reference Litt2019) (Fig. 2). We calculate daily melt over ice when the glacier surface is bare (using a TI of 
$6.5\,\mathrm{mm}^{\circ}\mathrm{C}^{-1}$ day−1, computed from values from Litt and others Reference Litt(2019) and melt from snow for pixels that are snow covered (using an index of 
$4.1\,\mathrm{mm}\,^{\circ}\mathrm{C}^{-1}$ day−1, following Braithwaite Reference Braithwaite(2008)). We scale the relative fraction of ice and snow melt per pixel using the relative snow cover data. While melt or surface runoff from rainfall from outside the model domain may also reach the model domain and eventually the glacier bed, we do not consider these inputs here. The TI model is shown to be more accurate for glaciers below 3500 m above sea level (a.s.l.) (Litt and others, Reference Litt2019), which is where most of Shishper’s tongue is located (Fig. 1). ERA5 data were downscaled from native 9 km to the model resolution (50 m) using a Kriging interpolation (Kusch and Davy, Reference Kusch and Davy2022). While some in situ climate data are available in the region, no station was operational in the vicinity of the glacier; using in situ data from an off-glacier station far away from the glacier would introduce its own set of uncertainties. Due to the relatively high temporal (1 day) and spatial (
$1\,\times\,1\,\mathrm{deg}$) resolutions, ERA5 data provide the best available estimate of meltwater inputs to the bed. SHAKTI is able to represent meltwater inputs as either point inputs or as distributed inputs; in this study, we apply the ERA5 melt estimate as a spatially distributed input over the bed, which is appropriate for heavily crevassed glaciers such as Shishper. The strong hydraulic coupling between surface and basal meltwater environments (Iken and Bindschadler, Reference Iken and Bindschadler1986; Zwally and others, Reference Zwally, Abdalati, Herring, Larsen, Saba and Steffen2002; Gulley and Benn, Reference Gulley and Benn2007; Shepherd and others, Reference Shepherd, Hubbard, Nienow, King, McMillan and Joughin2009; Miles and others, Reference Miles2017) has given us justification to make the assumption that all meltwater inputs to the bed (i.e. surface melt, rainwater, aquifer contributions) are instantaneous. Furthermore, the broken-up and crevassed nature of Shishper’s surface could allow for quicker delivery of meltwater to the bed.
Figure 2. (a) Englacial inputs to the transient subglacial hydrology model, averaged over the glacier, as calculated by ERA-5 Land and the temperature-indexed ablation model. (b) Average englacial input during the 2017 melt season (May through September). The red outline indicates the modeled domain.
3. Transient glacier hydrology simulations
3.1. Establishing winter base state
Before transient simulations can be run, the base winter state of the hydrological system must be established. To do this, we allow the drainage system to develop with zero external meltwater to the bed. During the winter, we assume that there is no surface or englacial melt, with geothermal flux and turbulent dissipation as the only sources of meltwater at the bed. Geothermal heat flux is set to 70 mW m−2, which is within previously measured values in the area (Shengbiao and Jiyang, Reference Shengbiao and Jiyang2000). Note that we have prescribed the sliding velocity to be zero, so there is no frictional heating or cavity opening from sliding over bumps. Because we exclude all contributions from tributary glaciers, a Neumann boundary condition of zero flux is applied to all lateral edges of the domain. A Dirichlet boundary condition is applied to the near-terminus domain boundary, with hydraulic head equal to the bed elevation (i.e. water pressure equal to atmospheric pressure at the outflow). A time step of 1 day is used for obtaining the final equilibrated state.
We define equilibrium by assessing the time rate of change of gap height, basal water flux, hydraulic head and effective pressure. We deem the model ‘equilibrated’ if there is no visible growth or decay in the minimum, maximum and spatial mean values of each of these parameters after 500 days; for example, mean effective pressure changes at a constant rate of 2e–7 % per year at the end of the winter equilibration. Once all output parameters reach equilibrium, after 600 days, there is formation of a primary drainage channel down the main trunk of the glacier (Fig. 3a). It is also worthwhile to note that the channel has formed in the absence of any surface water melt, indicating its potential to persist through the winter months just given a small amount of meltwater from geothermal flux and turbulent dissipation.
Figure 3. Basal flux across the modeled domain following (a) a ‘winter state’ equilibration spinup with no melt inputs to the system (b) a transient simulation through a full calendar year, including a summer melt season and return back to frozen winter conditions. (c) The difference between the two equilibrated states.
Beginning from the base winter state (Fig. 3a), we run a transient simulation of 1 year (January 1–December 31). Following this year, the model reaches a new stable winter state (Fig. 3b). The second stable state is largely similar to the first, but shows more efficient, concentrated drainage at a few areas, including the terminus and up-glacier at 4028 km N. Running additional melt seasons yields no additional changes in winter drainage patterns, indicating a new equilibrium. This perennial channel then forms the basis for the subglacial system during the melt season.
3.2. Seasonal evolution of subglacial hydrology
To understand how Shishper’s subglacial drainage network responds to seasonal changes in meltwater flux, we run transient simulations across a period of 5 years, 2017–21, using a timestep of 30 minutes. The transient input for these simulations is the temporally and spatially varying sum of ice melt, snow melt and liquid precipitation (Fig. 2), applied as distributed meltwater inputs to the subglacial system throughout our model domain of the main glacier trunk. Potential incoming melt inputs from tributary glaciers are not included.
Figure 4 illustrates changes in the configuration of the drainage system throughout 2017, which is representative of the pattern observed across all 5 years. We see a mostly closed system in winter (Fig. 4b), which transitions to a highly efficient, channelized system at the peak of the melt season (Fig. 4c). The basal flux mirrors the surface melt input trend, peaking around August (Fig. 2a), while hydraulic head and effective pressure stay mostly steady, apart from spikes at the beginning and end of the melt season. At the height of the melt season, the drainage system extends to the northernmost part of the domain, splitting into arborescent patterns characteristic of channelized drainage. By October 5, these channels then disappear, with the upper part of the system having completely shut down. Finally, the system returns to the winter state by late October (Fig. 4d and e).
Figure 4. (a) Model outputs for 2017, including hydraulic head, basal flux and effective pressure. Mean, min and max refer to the spatially averaged mean and the minimum and maximum values over the mesh. (b) Log 10 basal flux across the glacier at four times during the year: January 1 (winter), August 1 (peak melt), October 5 (drawdown of the drainage network at the end of the melt season) and October 18 (return to winter conditions).
The lower channel, which traverses the mid- to lower trunk clearly persists through every simulated winter in 2017–21, and can be seen in both images of the ‘closed’ state (Fig. 4a, d and e). We know that this channel appears during the winter equilibration, during which time the only water at the ice–bed interface comes from turbulent dissipation and geothermal heat flux. There is always a consistent stream of water, although small, that keeps the main channel open. Bhambri and others Reference Bhambri(2020) show that surface melt elevations move from 6400 m in peak summer to 3500 m at the end of winter (no surface melt is observed in December, January, or February) meaning that the bottom part of the glacier will always receive more melt, and is more likely to contain channels, than the top.
Figure 5 presents a closer look at the rapid decreases in effective pressure at the beginning of the melt season (May–July). As shown in Fig. 5a, coming from primarily distributed winter drainage with low transmissivity, the rise in hydraulic head due to the system’s inability to transport growing fluxes, and the rapid fall in head back to the equilibrium value, show that the system resolves spikes in water pressure by developing more efficient pathways (i.e. increasing transmissivity by opening new channels). The calculated spikes in hydraulic head in Figs. 4 and 5 are higher than realistic physical values, in which localized buildups of very high water pressure would more quickly be resolved through hydraulic jacking (local uplift where water pressure exceeds flotation) and/or fracturing of the overlying ice. Since these processes are not explicitly represented within SHAKTI, localized large water pressures may be resolved more slowly in the simulations, thus appearing as non-physical values. Setting a warmer ice temperature may also result in less extreme spikes (see Appendix C). Such a buildup of hydraulic head can be observed in Fig. 5b (June 5), where a large area of negative N (high water pressure) can be seen at the northern part of the domain. On June 10, this area has relaxed, and by June 30, the entire section has almost completely returned to the original state of effective pressure, around 2 MPa across the mesh. Figure 6 depicts the channel system that is established during and after these events, showing that an area of distributed, heavy flow around 4029 and 4030 km N quickly coalesces to a narrow and efficient channel in response to higher water pressures.
Figure 5. Top: (a) Melt input (blue) overlaid with spatially averaged hydraulic head (red) during 2017. Bottom: pressures across the mesh during (b) and after (c) the spike in hydraulic head in early June. Gray contours indicate effective pressure of 0 MPa (flotation, where ice overburden equals water pressure).
Figure 6. Basal fluxes surrounding an early-season spike (depicted in Fig. 5) show a transition at the upper trunk from distributed, sheetlike flow to efficient, channelized flow. The red arrows highlight the formation of a channel that occurs between June 5 and June 30.
So long as high fluxes continue, melt opening exceeds creep closure, keeping efficient channels open during the majority of the melt season. The drainage system is able to quickly shuttle large fluxes through, allowing it to return to a low-pressure state and draining the surrounding bed. Although velocities are not simulated here, it is inferred that sliding velocities would decrease due to a return to higher effective pressures in the summer. A velocity dataset developed by Beaud and others Reference Beaud, Aati, Delaney, Adhikari and Avouac(2022) at Shishper Glacier from 2013–19 shows that the glacier does indeed slow down significantly during summer months. In addition, increases in surface displacement further up the trunk of Shishper were observed by Bhambri and others Reference Bhambri(2020) during the early melt season (May to June) between 2013–16, suggesting that there could be decreased effective pressures at the northern part of the domain during this time. This agrees with our model results: near the terminus, the system remains perennially channelized, while the upper part sees an inefficient, distributed system during the early melt season that evolves to become more efficient over the summer.
As the system closes and the capacity of the drainage system falls, it regains its sensitivity to temporary increases in melt, as is seen in the early and late summer spikes in hydraulic head (Fig. 5) (Hart and others, Reference Hart, Young, Baurley, Robson and Martinez2022). This contraction happens as basal flux falls, allowing melt opening to fall and creep closure to dominate. The spikes are smaller than the ones at the beginning of the melt season because the system has not had much time to close yet, so it is more efficient than at the beginning of spring.
Overall, these findings corroborate the established understanding that there is a transition from a distributed to channelized drainage system and back during the course of the year (Fig. 4) (Schoof, Reference Schoof2010; Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013; Flowers, Reference Flowers2015). As long as high meltwater fluxes persist, melt opening exceeds creep closure, maintaining open channels throughout most of the melt season. As the system shuts down and the drainage capacity decreases toward the end of the melt season, it exhibits heightened sensitivity to melt increases, as evidenced in early and late summer hydraulic head spikes (Fig. 5).
4. Shishper surge phases between 2017 and 2019
4.1. Hydrological insights into surge dynamics
Comparing the modeled effective pressures with observed surge phases implies that incipient surge motion in November 2017 and subsequent slow acceleration through the winter 2017-18 do not show up as a clear hydrological signal in our simulations, suggesting that there could be a process or mechanism not accounted for by our model. However, significant hydraulic head spikes do correspond with rapid acceleration in June 2018, suggesting that elevated water pressures could play a role in escalating already-occurring ice motion (Kamb, Reference Kamb1987, Björnsson, Reference Björnsson1998).
Figure 7 overlays effective pressure simulated by SHAKTI on top of satellite-derived velocity observations from Beaud and others Reference Beaud, Aati, Delaney, Adhikari and Avouac(2022). Observations show a pre-surge acceleration begins in November 2017, but the model outputs indicate a decrease in effective pressure during this acceleration (Kamb, Reference Kamb1987; Björnsson, Reference Björnsson1998). At the beginning of June 2018, Bhambri and others Reference Bhambri(2020) describe a rapid but brief acceleration, which corresponds to the peak of 5.5 m d−1 described by Beaud and others Reference Beaud, Aati, Delaney, Adhikari and Avouac(2022) at the same time, coinciding with a series of modeled spikes in hydraulic head at the beginning of the 2018 melt season. As the drainage system enters its efficient summer state, the surge then enters a very slow ‘semi-quiescent’ period during which velocity is only slightly higher than normal summer velocities, lasting until September 2018 (Beaud and others, Reference Beaud, Aati, Delaney, Adhikari and Avouac2022). The glacier then accelerates again, reaching speeds of 2 m d−1 by November 2018 and 3.5 m d−1 in January 2019. Another surge peak occurs from late April to early May 2019 (Bhambri and others, Reference Bhambri2020). A small GLOF of the ice-dammed lake, which damaged the Karakoram Highway, follows from June 22–23, 2019 (Bhambri and others, Reference Bhambri2020; Beaud and others, Reference Beaud, Aati, Delaney, Adhikari and Avouac2022).
Figure 7. Glacier surface velocities from the dataset of Beaud and others Reference Beaud, Aati, Delaney, Adhikari and Avouac(2022) overlaid on model outputs of basal flux and effective pressure during the 2017–19 surge. The bright red line indicates a GLOF that occurred on June 22–23, 2019.
Our simulated spring and fall dips in effective pressure correspond with the observations by Beaud and others Reference Beaud, Aati, Delaney, Adhikari and Avouac(2022) of observations of spring and fall speedups at Shishper Glacier even during quiescent (non-surging) periods. The model results show larger and longer-duration effective pressure drops in spring compared to fall, aligning with observations of larger spring speedups. These spikes in hydraulic head appear more extreme than what may be expected in real life, in which localized buildups of very high water pressure would more quickly be resolved through hydraulic jacking and/or fracturing of the overlying ice. Overall, these findings support the hypotheses of observational studies suggesting that seasonal hydrology evolution is largely driving seasonal glacier motion trends in HMA (e.g. Nanni and others, Reference Nanni, Scherler, Ayoub, Millan, Herman and Avouac2023).
4.2. Limitations and future directions
Subglacial hydrology is likely only one of several factors that may drive surge behavior. To further disentangle the drivers of surge motion, it is necessary to consider and model additional processes such as frictional feedbacks due to sliding at the ice–bed interface, basal melting due to changes in the pressure melting point, till deformation, uplift and hydrofracture, dynamic advances and retreat of the terminus, and changes in ice thickness at the reservoir and receiving zone of the glacier (e.g. Flowers and others, Reference Flowers, Roux, Pimentel and Schoof2011; Liu and others, Reference Iken, Röthlisberger, Flotron and Haeberli2024). Furthermore, our model results hint that abrupt transitions from unstable, high water pressures to low (sub-flotation) pressures could play a role in slowdowns in surge motion. To quantify the role of subglacial hydrology in ice motion, a coupled model is necessary. Two-way coupling of SHAKTI with ice dynamics in ISSM has been implemented and applied recently to Helheim Glacier, Greenland (Sommers and others, Reference Sommers2024); implementing a similar coupled framework for this glacier could provide further insights into these complex interactions that are important for understanding the motion of surging glaciers. In addition, we have neglected the hydrological influence of upper-elevation tributary glaciers; although we do not expect large hydrological contributions from the tributaries due to their high elevation, the magnitude of hydrological flux from these tributaries may be non-trivial and ought to be considered in future work. Additional contributions from groundwater flow may also affect subglacial hydrology. Future studies should focus on integrating these processes into the model to better understand the interplay between subglacial hydrology, ice dynamics, ice-dammed lake floods, and surge behavior.
5. Summary and conclusions
Our study demonstrates that subglacial hydrology plays a crucial role in modulating glacier dynamics, particularly in surge-type glaciers like Shishper. The simulations show that at least one year’s melt cycle is required to bring the drainage system to a long-term equilibrium in which the subglacial drainage system returns to the same configuration every winter. This winter configuration features a primary channel in the lower trunk of the glacier, which remains year-round and serves as the basis for an arborescent, channelized drainage system that grows far up the glacier as the melt season peaks.
Our simulations demonstrate SHAKTI’s ability to represent the transition from an inefficient to efficient drainage pattern as melt flux rises and vice versa. These transitions are marked by large spikes in hydraulic head and corresponding dips in effective pressure, which support numerous previous observations of spring and fall speedups at Shishper and other mountain glaciers and strengthen existing hypotheses that seasonal glacier motion in HMA is largely driven by changes in subglacial hydrology. In addition, while subglacial hydrology is widely understood to be a crucial factor behind surging, the lack of a clear hydrological trigger for incipient surge motion and for the second surge peak highlights the complexity of surge dynamics and the need for further investigation into the interactions between subglacial hydrology, ice dynamics and other potential triggering mechanisms (Sevestre and Benn, Reference Sevestre and Benn2015, Benn and others, Reference Benn, Fowler, Hewitt and Sevestre2019).
This is the first time the SHAKTI model has been applied to a realistic mountain glacier. While our simulations here involve several simplifying assumptions to focus on the evolution of subglacial hydrology in isolation from velocity coupling, the model reproduces transitions between distributed and channelized drainage that could explain the timing of spring and fall speedups in the region. This provides a solid framework for future work to expand application of the model. These future studies should focus on analyzing the complex coupling between subglacial hydrology and glacier motion (Hoffman and Price, Reference Hoffman and Price2014, Sommers and others, Reference Sommers2024). Additionally, investigating the causal link between ice-dammed lake drainage and surge termination may provide valuable insights into the role of subglacial hydrology in modulating surge behavior (Björnsson, Reference Björnsson1998).
Data availability statement
ISSM (including SHAKTI) is freely available at issm.jpl.nasa.gov. Simulations were performed with ISSM version 4.24. Model output data and scripts are available in a Zenodo repository at https://doi.org/10.5281/zenodo.15644288 (Narayanan, Reference Narayanan2025). Plots in this paper make use of ColorBrewer colormaps developed by Stephen23 (2025).
Acknowledgements
We are grateful to Yoram Terleth, an anonymous reviewer, and Editors David Rounce and Hester Jiskoot for their feedback, which greatly improved the quality of the manuscript. We also thank Flavien Beaud for his helpful advice and velocity dataset, and Taylor Perron and Camilla Cattania for their insights. NN received support from the William Asbjornsen Albert Memorial Fellowship and the MIT Office of Graduate Education. NN and WC received support from the Heising-Simons Foundation HSF-2020-1909; ANS and CRM received support from HSF-2020-1911 and NASA grant 80NSSC24K1634.
Competing interests
The authors declare none.
Appendix A. Constants and parameter values
See Table A1.
Table A1. Constants and parameter values used in this study
Appendix B. Mesh resolution tests
We conducted a simple test of the finite element mesh resolution to ensure that the development of basal channels was not dependent on an arbitrary choice of mesh element size. We ran winter equilibrations with triangular mesh sizes of 10, 20, 40, 50, 100, 200 and 250 m. Each was run with 6-hour timesteps for 300 days and was initialized with the same initial conditions. In Fig. B1, we show gap heights at the end of each of these winter equilibrations.
Figure B1. Gap height (m) across the domain, shown for mesh resolutions of 10, 20, 40, 50, 100, 150 and 200 m.
Areas of high gap height show the location of channels and subglacial lakes. The location of these channels is largely invariant with mesh resolution, suggesting that channel locations exhibit a higher dependence on topography than on mesh resolution.
Appendix C. Sensitivity to creep parameter A
Since the temperature at the base of Shishper is unknown, we conduct a brief sensitivity test to assess how results may change given different basal ice rheologies. In most of the simulations featured in this study, we use a value for Glen’s flow law parameter A corresponding to a temperature of 
$-5^{\circ}\mathrm{C}$ (see Appendix A). Here we re-run 2017 with identical forcings but using an A value corresponding to temperate ice (
$0^{\circ}\mathrm{C}$). Figure C1 shows spatially averaged means corresponding to 
$-5^{\circ}\mathrm{C}$ ice and 
$0^{\circ}\mathrm{C}$ ice and the difference between the two.
Figure C1 demonstrates that the timing of transitions between efficient and inefficient drainage remains the same regardless of ice softness. The amplitude of basal flux and hydraulic head remain largely the same as well; however, the amplitude of effective pressure has slightly more significant differences between the two simulations. Notably, warmer ice corresponds with less extreme spikes in effective pressure; the largest differences between the simulations occur during the early melt season spikes.
Figure C1. Mean values for (a) effective pressure, (b) basal flux and (c) hydraulic head. The first simulation (Nmean1, qmean1 and hmean1) was modeled with A corresponding to 
$-5^{\circ}\mathrm{C}$ while the second simulation (Nmean2, qmean2 and hmean2) used A corresponding to temperate ice (
$0^{\circ}\mathrm{C}$). ‘Diff’ refers to the difference between the two simulations.
Open access
