1. Introduction
Conditions at the beds of ice sheets and glaciers strongly influence ice dynamics through the effect of lubrication and enhanced sliding where basal water pressure is high. Unfortunately, direct observations are difficult to obtain beneath hundreds to thousands of meters of ice. Techniques such as drilling boreholes to the bed (Iken and others, Reference Iken, Echelmeyer, Harrison and Funk1993; Murray and Clarke, Reference Murray and Clarke1995; Harper and others, Reference Harper, Humphrey, Pfeffer, Fudge and O'Neel2005; Fudge and others, Reference Fudge, Humphrey, Harper and Pfeffer2008; Andrews and others, Reference Andrews2014; Ryser and others, Reference Ryser2014), using radar sounding to infer the presence of liquid basal water (Oswald and Gogineni, Reference Oswald and Gogineni2008; Chu and others, Reference Chu2016; Jordan and others, Reference Jordan1999; Oswald and others, Reference Oswald, Rezvanbehbahani and Stearns2018), and dye-tracing tests (Nienow and others, Reference Nienow, Sharp and Willis1998; Cowton and others, Reference Cowton2013) are helpful means to gain a view into basal conditions, but do not provide a complete description of the spatially and temporally heterogeneous bed environment.
Several numerical models have been developed to simulate the flow and pressure of water beneath glaciers and ice sheets (e.g., Flowers, Reference Flowers2015; de Fleurian and others, Reference de Fleurian2018) and have successfully reproduced melt-season drainage system evolution. However, challenges remain in these efforts, and subglacial hydrology model development remains an active area of research. One persistent issue is that many models rely on unconstrained parameters, for example, prescribing a typical height and spacing of asperities at the bed, or specifying hydraulic conductivity of drainage system components. Subglacial hydrology models are sensitive to these uncertain parameters, with small changes leading to substantial differences in simulated basal water pressures and drainage regimes (Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013; Banwell and others, Reference Banwell, Hewitt, Willis and Arnold2016).
Another challenge is that models have had difficulty reproducing widespread areas of high winter water pressures (Flowers, Reference Flowers2015) that have been observed in Greenland borehole arrays (Harper and others, Reference Harper, Humphrey, Pfeffer, Fudge and O'Neel2005; Ryser and others, Reference Ryser2014). Recent work highlighting the importance of hydraulically disconnected regions of the bed and incorporating a representation of these isolated areas into drainage models has helped to explain this phenomenon (Hoffman and others, Reference Hoffman2016; Rada and Schoof, Reference Rada and Schoof2018; Rada Giacaman and Schoof, Reference Rada Giacaman and Schoof2023), linking disconnected regions with other drainage components that represent flow through inefficient sheet-like configurations and efficient channels. Different drainage ‘modes’ are treated with disparate equations to represent distinct physical processes (for example, channels open by melting while the sheet-like system opens by sliding over asperities in the bed). However, distinguishing between drainage modes by applying separate equations to represent different components imposes an artificial distinction and may fail to fully capture the continuous spectrum of spatio-temporally evolving drainage behavior that exists in reality. In addition, much of the previous work modeling subglacial hydrology in Greenland focuses on land-terminating glaciers in west Greenland. Although recent modeling efforts have been applied to Greenland tidewater glaciers such as Store Glacier (Cook and others, Reference Cook, Christoffersen, Todd, Slater and Chauché2020, Reference Cook, Christoffersen and Todd2022), Petermann Glacier (Ehrenfeucht and others, Reference Ehrenfeucht, Morlighem, Rignot, Dow and Mouginot2023), and Helheim Glacier (Stevens and others, Reference Stevens2022), not as much attention has been directed at tidewater glaciers in this context, despite their outsize influence on mass loss from the ice sheet.
In this paper, we describe a reduced form of the Subglacial Hydrology And Kinetic, Transient Interactions (SHAKTI) model to address the above-mentioned issues. First presented by Sommers and others (Reference Sommers, Rajaram and Morlighem2018), SHAKTI takes a continuum approach without explicitly distinguishing between different drainage components, yet does represent behavior corresponding to different ‘modes’ of drainage, primarily facilitated through flow regime transitions with a single set of equations. Here, we eliminate some terms that rely on unconstrained parameters or uncertain physics, and which are not needed in order to obtain a minimal model. After summarizing the equations governing the evolution of the subglacial hydrology system in SHAKTI, we apply the model to Helheim Glacier in east Greenland under winter conditions to demonstrate its capabilities in a real glacial setting and attempt to reconcile the outstanding problem of simulating high winter water pressures with a continuum model. We also examine the role of basal melt from frictional heat and the influence of different styles of calculating basal shear stress.
2. Methods
2.1 Summary of equations
SHAKTI uses a single set of equations to calculate hydraulic head, effective pressure, basal water flux, and geometry of the subglacial drainage system. In contrast to other subglacial hydrology models, SHAKTI allows for natural transitions between laminar and turbulent flow, allowing distinct flow regimes to coexist in different regions of the model domain with spatially and temporally variable transmissivity, giving rise to a spectrum of inefficient and efficient drainage configurations. SHAKTI includes heat generated by energy dissipation within the water flow and opening of the drainage system by melt across the entire domain, unlike models that treat ‘inefficient’ sheet-like and ‘efficient’ channel-like components of the drainage system with different equations. SHAKTI's unified approach leads to the emergence and disappearance of flexible drainage configurations over time, conserving mass and energy within the system. In what follows, we summarize the original SHAKTI model equations, along with key modifications to eliminate terms that depend on unconstrained parameters. Whereas in the original SHAKTI formulation, we included terms to facilitate direct comparison to other models (de Fleurian and others, Reference de Fleurian2018), here we examine the model capabilities in the absence of these unconstrained terms. Tables 1 and 2 serve as convenient references for the variables, constants, and parameters used in the equations. For a complete description of the original SHAKTI model, we refer readers to Sommers and others (Reference Sommers, Rajaram and Morlighem2018).
Variables
Constants and parameters used in this study
SHAKTI is composed of partial differential equations that describe conservation of ice and water mass, drainage configuration, water flux, and internal melt generation. The water balance equation is written as
where b is the gap height between the ice and bed, b e refers to a volume of water stored englacially per unit area of the bed, q is gap-integrated water flux through the subglacial system, $\dot m$
is the melt rate expressed as a mass flux (units of kg m−2 s−1), ρ w is density of water, and i e→b is the rate of meltwater input from the englacial system to the bed, which can be specified and handled by the model as a combination of distributed input (units of m s−1) and point inputs to represent moulins or crevasses (units of m3 s−1).
For our current purposes of simulating a steady-state minimal winter base state, the term ${\partial b_e\over \partial t}$
is zero. This implies that all simulated water is at the bed and we do not attempt to approximate delayed release from englacial storage, although water held in void spaces within the ice could play an important role in subglacial water flow even in the absence of surface melt (Schoof and others, Reference Schoof, Rada, Wilson, Flowers and Haseloff2014). For winter conditions, with no englacial storage and no external meltwater inputs, the modified water balance equation is then given by
The geometry of the drainage system is represented by the average gap height b over a discrete area of the bed, which evolves through time dynamically. In the original equations, gap height increases as a result of both melt and by sliding over bumps in the bed, and decreases due to creep deformation. This can be expressed as change in gap height over time as
where ρ i is the density of ice, β is a dimensionless coefficient that dictates opening of the subglacial gap by sliding over bumps in the bed, ub is the ice sliding velocity, A is the flow law parameter, p i = ρ i g H is ice overburden pressure in which g is gravitational acceleration and H is ice thickness, p w = ρ w g (h − z b) is subglacial water pressure in which h is hydraulic head and z b is bed elevation above sea level, and n is Glen's flow law exponent.
The coefficient β that governs opening by sliding over bumps in the bed depends on prescribing an uncertain ‘typical bed-bump height’ (b r) and ‘typical bed-bump spacing’ (l r), following Werder and others (Reference Werder, Hewitt, Schoof and Flowers2013). This method of cavity opening was introduced by Hewitt (Reference Hewitt2011) to parameterize an opening mechanism in distributed drainage other than melt, based on the description of linked cavities of Kamb (Reference Kamb1987). Schoof and others (Reference Schoof, Hewitt and Werder2012) found that a system of linked cavities that opened by melt was unstable (assuming turbulent flow); therefore, the opening-by-sliding mechanism has been widely adopted for the evolution of inefficient drainage systems. However, to prevent numerical instability in other models that use this type of opening mechanism, sliding velocity must usually be capped. For example, Poinar and others (Reference Poinar, Dow and Andrews2019) applied the Glacier Drainage System model (GlaDS) (Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013) with b r = 0.08 m and l r = 2 m to an idealized Helheim Glacier-like domain, and had to cap sliding speed at 800 m yr−1 to achieve model stability. This is an order of magnitude less than observed velocity on Helheim Glacier. Stevens and others (Reference Stevens2022) used a similar model to simulate subglacial hydrology at Helheim and chose a bed bump height of 0.1 m and spacing of 10 m without imposing a limit on the velocity, but did not address sensitivity of the resulting subglacial drainage system to these choices of bed roughness geometry. Applying GlaDS to Sermeq Kujalleq (Store Glacier), a tidewater glacier in west Greenland, Cook and others (Reference Cook, Christoffersen, Todd, Slater and Chauché2020, Reference Cook, Christoffersen and Todd2022) selected higher bed bump and spacing values, b r = 1 m and l r = 100 m. In fast-flowing glaciers, we find that including opening by sliding over bumps in SHAKTI effectively smooths out and inhibits any channelized structure, leading to an unrealistic near-uniform gap height equal to the specified bed bump height over large regions. This opening mechanism was included in the original SHAKTI equations largely for comparison with other similar models, but is not required for model stability due to the transitional flux formulation that allows for changes between laminar and turbulent flow regimes, discussed below in Equation (5) and by Sommers and others (Reference Sommers, Rajaram and Morlighem2018).
Given how the opening-by-sliding parameterization depends on arbitrarily prescribed bed-bump height and spacing that in reality are heterogeneous, we cannot be confident that the commonly used formulation accurately represents the increase in average b – especially in the case of fast-flowing glaciers, which we expect to be underlain by till. We eliminate the ‘opening by sliding over bumps at the bed’ term β|ub| and allow the drainage geometry to open only due to melt, everywhere in the domain, which behaves stably with our transitional flux formulation (Eqn. (5) below). The minimum gap height allowed is 10−3 m, representing a transition to premelted films (Wettlaufer and Worster, Reference Wettlaufer and Worster2006; Rempel and others, Reference Rempel, Meyer and Riverman2022). We write our modified basal gap dynamics equation as
The momentum equation that describes the water flux is
where ν is kinematic viscosity of water, ω is a parameter related to a friction factor that controls the transition from laminar to turbulent flow, and Re = |q|/ν is the local Reynolds number. When ωRe ≪ 1, Eqn. (5) behaves like laminar flow, with q proportional to the magnitude of the hydraulic gradient $\vert \nabla h\vert$
. When ωRe ≫ 1, by inserting the definition Re and rearranging to solve for q, we see that q is proportional to $\sqrt {\vert \nabla h\vert }$
, corresponding to completely turbulent flow. This flux formulation is based on equations developed in the context of flow in rock fractures (Zimmerman and others, Reference Zimmerman, Al-Yaarubi, Pain and Grattoni2004; Chaudhuri and others, Reference Chaudhuri, Rajaram and Viswanathan2013; Rajaram and others, Reference Rajaram, Cheung and Chaudhuri2009). The general Eqn. (5) also allows for intermediate transitional regimes. This flux formulation or ‘flow law’ is the key difference of the original SHAKTI model compared to other subglacial hydrology models, and plays an important role in maintaining stability while allowing for channelization to occur. Forcing q to be always laminar or always turbulent (by changing the value of ω) results in a model instability in some situations, but these scenarios behave well when allowing for the flexible flow regime transition around Re = 103 (ω = 0.001), accordingly generating spatially variable Re distributions (Sommers and others, Reference Sommers, Rajaram and Morlighem2018).
In contrast to models that rely on a prescribed hydraulic conductivity, this flux formulation incorporates a spatio-temporally variable hydraulic transmissivity K, given by
Basal meltwater generation is represented in SHAKTI through an energy balance at the bed, assuming ice and water are both always at the pressure melting point, i.e.
where L is the latent heat of fusion of water, G is geothermal heat flux, τ b is the basal stress, c t is the change of pressure melting point temperature with pressure, and c w is the heat capacity of water. This energy equation includes melt due to heat contributed by geothermal sources, frictional heat from sliding over the bed, turbulent dissipation, and adjustments for changes in the pressure melting point due to changes in water pressure. Geothermal flux varies spatially, and the value we use here (see Table 2) corresponds to the lower-end estimate of values typical for Greenland (Martos and others, Reference Martos2018). Note that Equation (9) in Sommers and others (Reference Sommers, Rajaram and Morlighem2018) should have a positive sign for the last term as written here in Equation (7), which accounts for changes in the sensible heat due to changes in the pressure-melting-point temperature with changes in water pressure. This term is generally a modest heat sink for flat beds, reducing the melt rate as in Röthlisberger (Reference Röthlisberger1972), but can also contribute to enhanced melt with steep slopes or supercooling when flowing uphill in bed overdeepenings (see, e.g. Creyts and Clarke, Reference Creyts and Clarke2010). In regions of the bed where ice is separated from the bed by a thin film, however, the difference between ice normal stress and water pressure will be balanced by premelting, where the relationship between pressure and temperature does not follow the Clapeyron slope (Rempel and Meyer, Reference Rempel and Meyer2019). Co-located observations of pressure and temperature within subglacial boreholes show that the Clapeyron relationship does not always hold (Huss and others, Reference Huss, Bauder, Werder, Funk and Hock2007). We eliminate this term for now in the reduced SHAKTI model, and leave further evaluation of the effects of pressure melting on subglacial hydrology to future work.
Basal shear stress as implemented here depends on a drag coefficient C, effective pressure N, and sliding velocity ub, τ b = C 2N|ub| (Budd and others, Reference Budd, Keage and Blundy1979). We use observed surface velocity as a proxy for basal sliding velocity, a reasonable assumption in fast-flowing glaciers. We also examine results using other formulations for τ b to explore the impact of how frictional heat from sliding is represented.
The melt rate at the bed considered here in the reduced form of SHAKTI is a result of geothermal flux, frictional heat from sliding, and heat generated by mechanical energy dissipation in the subglacial system:
We combine equations (2), (4), (5), (6), and (8) to form an elliptic equation in terms of hydraulic head:
We solve Eqn. (9) for the head distribution using a Picard iteration to handle the nonlinear dependence of the terms on the right-hand side of the equation, then we solve Eqn. (4) explicitly to evolve the gap height b. No numerical limits are imposed on head (i.e., water pressure is free to exceed overburden pressure or to become negative over the course of a simulation).
SHAKTI is built into the Ice-sheet and Sea-level System Model (ISSM; Larour and others, Reference Larour, Seroussi, Morlighem and Rignot2012) using the finite element method in a parallelized computational framework. In addition to the elimination of terms that rely on uncertain parameters, we modify how the term that describes creep closure in Equation (9) is handled within the Picard iteration. This change helps the Picard iteration converge instead of oscillating, a problem that arises under thick ice with low meltwater input. In previous work with SHAKTI (de Fleurian and others, Reference de Fleurian2018; Sommers, Reference Sommers2018), this oscillation obstacle was handled using under-relaxation. With a Newton linearization weighting, inspired by Gagliardini and Werder (Reference Gagliardini and Werder2018) in their implementation of a similar subglacial hydrology model, GlaDS (Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013) in a different ice-sheet model, Elmer/Ice (Gagliardini and others, Reference Gagliardini2013), this change to the creep term numerics facilitates convergence of the iterative process to find h. This is a key practical improvement for the application of SHAKTI to glacial environments with thick ice and low water inputs, common in Greenland during the winter.
2.2 Model domain
We explore the winter base-state subglacial drainage of Helheim Glacier in east Greenland (Fig. 1). Helheim is a large, fast-flowing glacier that terminates in Sermilik Fjord with two main branches that flow through deeply incised canyons. We use bed and ice surface elevation based on BedMachine (Morlighem and others, Reference Morlighem2017). Our model domain includes the two main branches of Helheim Glacier and extends inland to approximately 1900 m surface elevation. Other recent work simulating subglacial hydrology at Helheim by Stevens and others (Reference Stevens2022) explored a smaller domain focused near the terminus region, including the smaller tributary to the north, and found the small north tributary to be relatively insignificant in winter. Our domain is discretized into an unstructured triangular mesh refined based on ice velocity, with 12472 finite elements and 6,371 vertices (Fig. S1). The area of each element ranges from $1.00\times 10^4 - 4.39\times 10^6\, {\rm m}^2$
($0.01 - 4.39$
km2), with a mean area of 1.65 × 105 m2 (0.165 km2).
(a) Location of Helheim Glacier on the Greenland ice sheet (inset), velocity (Joughin and others, Reference Joughin, Smith and Howat2018) over model domain used in SHAKTI simulations, overlaid on 2010 MODIS mosaic (Haran and others, Reference Haran, Bohlander, Scambos, Painter and Fahnestock2018), (b) bed topography in model domain relative to sea level (Morlighem and others, Reference Morlighem2017, Reference Morlighem2021), (c) surface elevation in model domain relative to sea level.
2.3 Boundary conditions
We set the boundary condition for hydraulic head at the glacier terminus as a Dirichlet condition, based on the ideas that the terminus is grounded and that pressure at the subglacial outflow should be equal to hydrostatic pressure from the overlying fjord water (depth varying across the front). Setting the subglacial water pressure equal to the pressure in the fjord at the subglacial outflow gives
where ρ f is the density of the fjord water and d is the water depth (d = −z b is a positive quantity, where z b is bed elevation relative to sea level, i.e. a negative quantity at the terminus). Rearranging, we solve for head:
where ρ f/ρ w > 1. If the water at the subglacial outflow is assumed to be well-mixed with fresh water from melting at the glacier front, then ρ w ≈ ρ f, and therefore h ≈ 0 at the outflow. In other recent simulations of Helheim subglacial hydrology, Stevens and others (Reference Stevens2022) assumed a floating terminus.
We prescribe Neumann conditions on the upstream and lateral boundaries of the model domain, with $\nabla h = {\bf 0}$
at these boundaries.
2.4 Experiment description
Using the reduced form of the SHAKTI model as described above, we conduct simulations to produce a representation of the minimal winter base-state hydrological system of Helheim Glacier. In winter, no surface meltwater is produced, and meltwater inputs that reach the bed from the surface are presumed to be essentially zero over most of the Greenland Ice Sheet. Accordingly, we assume the absence of discrete features (e.g. moulins) or delayed drainage from features such as firn aquifers, in contrast to Poinar and others (Reference Poinar2017, Reference Poinar, Dow and Andrews2019). This does not mean, however, that there is no water at the bed in winter, as water is generated at the glacier base (Eqn. (8)). To generate an estimate of the winter base state of the subglacial hydrological system of Helheim Glacier, we run a ‘spin-up’ SHAKTI simulation with zero meltwater input from the englacial system to the bed (i e→b = 0) with all water at the system produced by basal melt as calculated by Eqn. (8). The simulation is run for a period of 1 year, with a time step of 30 minutes. Figure S2 presents time series of various output quantities to demonstrate convergence of the system within this simulation duration.
In addition to this winter spin-up simulation, which we will refer to as main, we perform additional simulations using the same model domain, time step, and duration, to test and demonstrate the influence of frictional heat. We explore three different approaches of calculating the basal shear stress τ b (involved in the frictional heating contribution to basal melt): 1) basal shear stress equal to the driving stress, ${\tau }_b = \rho _igH\vert \nabla z_s\vert$
, 2) basal shear stress equal to yield stress, depending on evolving effective pressure, τ b = μN, where μ is the till friction coefficient (Iverson, Reference Iverson2010), and 3) zero basal shear stress, τ b = 0, which effectively removes the influence of frictional heat from sliding.
Table 3 summarizes the simulations presented in this paper. The focus of this study is winter hydrology, and transient simulations will be the focus of a separate publication to maintain this focus. However, example simulation results with time-dependent forcing are shown in the supplementary material (Figs. S6 and S7) to illustrate the robustness of the reduced SHAKTI model with realistic topography under seasonal conditions.
Summary of simulations
3. Results
3.1 Winter base state
In Figure 2, we present the spatial distribution of the subglacial drainage system at the completion of the main simulation base-state spin-up. In this winter base state with zero external meltwater input, large portions of the bed exhibit high water pressure as demonstrated by a fraction of overburden, p w/p i, that is close to one (Fig. 2a). Distinct preferential drainage pathways emerge with larger gap height (Fig. 2b) and higher Reynolds number and water flux (Re, Fig. 2c) forming river-like structures. The major river-like structures coincide with the locations of deeply incised bedrock channels. Effective pressure is highly variable across the domain, lowest near the terminus and lower in the main drainage pathways than in the surrounding bed (Fig. 2d). Regions of turbulent flow (Re > 103) as well as regions of laminar flow (Re < 103) coexist (Fig. 2c), with clearly higher Re in the main pathways and in smaller arborescent tributaries that feed them. Most of the bed away from the main pathways has very low Re (i.e., very low flux, since $\vert q\vert = Re\ \nu$
), with regions that appear to be poorly connected. A clear primary outflow structure emerges at the terminus, located slightly south of the center line. This preferential discharge location agrees well with the location of observed summertime subglacial plumes upwelling at Helheim (Everett and others, Reference Everett, Murray, Selmes, Holland and Reeve2021; Melton and others, Reference Melton2022) and with the consistent outflow point as simulated by Stevens and others (Reference Stevens2022).
Winter base state of subglacial hydrological system with zero external meltwater input (main simulation): (a) water pressure as fraction of overburden, p w/p i, (b) gap height (shown in log10 scale for detail), (d) Reynolds number and basal water flux (shown in log10 scale for detail), (d) effective pressure, (e) transmissivity (shown in log10 scale for detail), (f) basal melt rate (shown in log10 scale for detail).
3.2 Spatially variable transmissivity
In recent years, the community has highlighted the importance of hydraulically isolated and weakly connected regions of the bed in subglacial hydrology, particularly for maintaining observed high winter water pressures (Andrews and others, Reference Andrews2014; Hoffman and others, Reference Hoffman2016; Rada and Schoof, Reference Rada and Schoof2018; Mejía and others, Reference Mejía2021; Rada Giacaman and Schoof, Reference Rada Giacaman and Schoof2023). An advantage of SHAKTI is that the flux formulation (Equation (6)) incorporates spatially and temporally variable hydraulic transmissivity, rather than requiring a prescribed value for hydraulic conductivity as in other models. As shown in Figure 2e, we find highly heterogeneous transmissivity values, varying over several orders of magnitude within the Helheim domain. Simulated transmissivity is highest near the terminus and along the main winter drainage pathways, and is particularly low along the divide between the two main branches of fast ice flow, as well as through the center of each of these ice streams. The low-transmissivity regions in the center of the main ice flow branches coincide with topographical ridges in the bed, particularly in the northern branch. These low-transmissivity areas represent regions with little connectivity and water flow through them, or in other words are interpreted as poorly connected. In groundwater aquifers, transmissivity of K < 5 m2 d−1 is considered to be negligible, 5 < K < 50 m2 d−1 is weak, 50 < K < 500 m2 d−1 is moderate, and K > 500 m2 d−1 is high (De Wiest, Reference De Wiest1965). Considering the weak transmissivity threshold of K ≤ 50 m2 d−1, 71% of the bed by area in our model domain is interpreted to be poorly connected in the winter.
The full range of hydraulic transmissivity shown in Figure 2e spans a range of 4.3 × 10−4 to 5.9 × 103 m2 s−1. To put our simulated transmissivity values into context of some other subglacial hydrology studies, Stevens and others (Reference Stevens2022) simulated partially overlapping range of transmissivity of 9.8 × 10−5 to 9.8 m2 s−1 (based on reported transmissivity values for one point of T = 1 × 10−8 to 1 × 10−3 Pa−1 s−1 with transmissivity defined in terms of pressure gradient, converted to standard transmissivity units of m2 s−1 as Tρ w g for comparison to our transmissivity dependent on hydraulic head gradient). Inferring transmissivity from GPS measurements in response to lake drainage events, assuming Darcy flow through porous media at the bed, Lai and others (Reference Lai2021) estimated transmissivity of 4.4 × 10−3 to 1.2 m2 s−1 (based on reported kh 0 values of 0.8 to 215 mm3, converted to units of m2 s−1 as kh 0g/ν). In modeling a land-terminating region in west Greenland, Dow and others (Reference Dow2015) calculated transmissivity of 2.5 × 10−3 to 0.2 m2 s−1, with a ‘baseline’ value of 1.5 × 10−2 m2 s−1, based on different values for sheet ‘conductivity’ (units of m s−1) and critical sheet thickness (units of m).
Our simulated transmissivity spans a much wider range than what has been previously inferred. To examine the spatial distribution of transmissivity and relation to water pressure, Figure 3 presents transmissivity as a function of water pressure (as fraction of overburden), grouped by region designated by surface elevation. The broadest range of transmissivity emerges in the low-elevation region, with z s ≤ 900 m. The intermediate elevation class, between 900 and 1500 m surface elevation, inhabits a middle range of transmissivity, while the highest elevation region in the interior (z s ≥ 1500 m) is broadly characterized by high fraction of overburden and very low transmissivity.
Winter hydraulic transmissivity (K) as a function of water pressure (as a fraction of overburden, p w/p i) for different regions of the model domain based on surface elevation.
3.3 Basal melt
Basal melting has increasingly been acknowledged as an important consideration for ice dynamics (Karlsson and others, Reference Karlsson2021; Young and others, Reference Young2022). As described in Equation (8), basal meltwater is produced at the bed through geothermal flux, frictional heat from sliding of the ice over the bed, and turbulent dissipation, in which mechanical energy is converted to thermal energy in the water flow. The mean melt rate over the entire domain is 0.097 m yr−1 and the total melt rate over the domain is 6.3 m3 s−1. Previous work exploring subglacial hydrology at Helheim or on an idealized Helheim-like domain using other subglacial hydrology models invoke a prescribed uniform background basal melt rate to account for geothermal and frictional heat (Poinar and others, Reference Poinar, Dow and Andrews2019; Stevens and others, Reference Stevens2022). In Figure 4, we present the fraction of basal melt rate due to geothermal flux (Fig. 4a), frictional heat (Fig. 4b), and dissipation (Fig. 4c). Geothermal flux is applied in our simulation with a uniform value of 0.05 W m−2, but in reality this varies spatially and may vary by up to a factor of two in narrow incised canyons (Colgan and others, Reference Colgan2021; Willcocks and others, Reference Willcocks, Hasterok and Jennings2021). Frictional heat dominates over most of the domain, yielding the highest melt rates along the steep topographic walls of the deeply incised bed canyons. Dissipation is an important source of basal melt in the primary outflow to the terminus, even in winter.
(a) Fraction of basal melt rate due to geothermal flux, (b) fraction of basal melt rate due to frictional heat from sliding, (c) fraction of basal melt rate due to dissipation.
In the results presented above, the basal stress τ b = C 2N|ub|, where τ b evolves transiently with N and basal velocities are held constant in time. The spatially variable drag coefficient, C, shown in Figure 5a, is obtained through inverse modeling using ISSM, assuming effective pressure N = ρ igH + ρ wgz b, the standard practice for basal drag inversion in ISSM (although the distribution of N obviously evolves in our SHAKTI simulations). The inversion optimizes C in order for the ice flow model to reproduce observed surface velocities through the ice stress balance. As we might expect intuitively, the resulting drag coefficient from inversion is high in areas of slow-moving ice (high friction) and low where the ice is sliding rapidly (high slip rates, i.e. in the main glacier branches and near the terminus). Basal shear stress at the end of the main simulation (Fig. 5b) varies over many orders of magnitude, with the highest stress along the walls of the main branches and throughout the interior, with lowest stress near the terminus and particularly along the primary discharge pathway.
(a) Drag coefficient used in main simulation basal stress calculation, (b) winter basal shear stress, τ b = C 2N|ub|.
Now we consider the three winter simulations with different approaches to calculate the basal shear stress τ b: 1) basal shear stress equal to the driving stress, ${\tau }_b = \rho _igH\vert \nabla z_s\vert$
(drivingstress), 2) basal shear stress equal to yield stress, depending on evolving effective pressure, τ b = 0.3N, where 0.3 is the till friction coefficient (yieldstress), and 3) zero basal shear stress, τ b = 0, which effectively removes the influence of frictional heat from sliding (nofrictionheat).
The resulting winter base-state hydrological system differs substantially between these three cases and our original main winter simulation. Overall, effective pressure is higher over most of the domain (corresponding to lower water pressure) when frictional heat is included than for the case with τ b = 0 (Fig. 6a). Drivingstress and yieldstress produce higher basal stress than our main simulation (Fig. 6b), and correspondingly lead to higher melt rates. With basal shear stress equal to driving stress or yield stress, frictional heat from sliding is the vastly dominant source of basal melt rate over most of the domain and leads to much higher flux through the main branches of the glacier, blurring the distinct drainage pathways that emerge in our main simulation and in the nofrictionheat simulation, as illustrated by the distribution of basal water flux in each simulation (Fig. 7). In our main simulation, water flux and basal melt rates are high along the walls of the deeply incised bed troughs (Figs. 2f and 7a). Better defined flow paths are visible than in the drivingstress and yieldstress simulations, but the drainage regime is distinctly more developed with a higher flux of water compared to the nofrictionheat case.
Histograms comparing different approaches for basal shear stress τ b: (a) resulting effective pressure distribution, (b) basal shear stress.
Winter basal water flux (shown in log10 scale for detail) resulting from different approaches for basal shear stress τ b: (a) drag coefficient from ISSM inversion, (b) no frictional heat, (c) driving stress, (d) yield stress.
We find that total basal melt rates and subglacial discharge at the terminus vary substantially depending on the frictional heat calculation. The total basal melt rate and outflow at the terminus with frictionheat are lower than in our main simulation, and several times higher with yieldstress and drivingstress, as summarized in Table 4.
Subglacial discharge at terminus and total basal melt rate over model domain
4. Discussion
4.1 High water pressure in winter
While models have achieved good qualitative behavior for melt-season evolution (Hewitt, Reference Hewitt2013; Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013), a challenge of subglacial hydrology modeling has been to reproduce high water pressures in winter conditions to agree with borehole measurements (Flowers, Reference Flowers2015). Note that most borehole measurements of high winter water pressure have been taken in land-terminating portions of west Greenland, with no borehole measurements to the bed at Helheim Glacier. Disconnected, weakly connected, or isolated regions of the bed have been shown to be necessary for maintaining high winter water pressure in other subglacial hydrology models (Hoffman and others, Reference Hoffman2016; Rada and Schoof, Reference Rada and Schoof2018). Our winter SHAKTI simulation of Helheim successfully produces widespread high water pressure (Figs. 2a, S3) that corresponds with low transmissivity in the interior (Fig. 2e), suggesting that the transitional flux formulation of SHAKTI (Equation (5)) enables representation of winter high-pressure regions with a continuum approach.
Recently, Stevens and others (Reference Stevens2022) simulated subglacial hydrology at Helheim with a different model (Hewitt, Reference Hewitt2013; Stevens and others, Reference Stevens, Hewitt, Das and Behn2018). While their winter simulation produced high water pressure (low effective pressure) in the main trunk of the glacier near the terminus, the interior portions of their domain extending up the two main branches were simulated to exist at effective pressure of N = 6 MPa or more during the winter, which is significantly higher than the maximum effective pressure of N = 3.2 MPa in our main simulation. Key differences in our approach compared to Stevens and others (Reference Stevens2022) that likely play a role in this difference are the transitional flux (Equation (5)) used in SHAKTI, and spatially variable basal stress.
4.2 Role of basal topography
The predictions for winter subglacial water pressure in Figures 2a and S3 show variability across a range of length scales. Hydropotential flow routing (Shreve, Reference Shreve1972) is commonly used to predict possible subglacial paths. This is typically done by assuming the water pressure to be equal to overburden everywhere (p w/p i = 1.0), or some other uniform fraction of overburden. However, as shown in Figure 2a, we find that the fraction of overburden varies considerably, spanning the entire range from 0 to 1 over the whole domain, and spanning from ~0.7 − 1 in regions of faster ice flow. Assuming a uniform fraction of overburden may thus yield an inaccurate picture of subglacial drainage (Wright and others, Reference Wright, Siegert, Le Brocq and Gore2008). In Figure 8, we compare the difference in head distributions from a uniform fraction of overburden and our winter SHAKTI simulation. An assumption of water pressure equal to 100% overburden pressure (p w = p i, i.e. a uniform fraction of overburden p w/p i = 1.0) overpredicts head in the vast majority of the domain, except very near the terminus where it agrees with the results of our SHAKTI simulation (Fig. 8a). Assuming water pressure equal to 80% overburden pressure (p w = 0.8p i, i.e. a uniform fraction of overburden p w/p i = 0.8) agrees better, but still overestimates head in the more stagnant portions of the domain lateral to the main drainage pathways and underestimates head near the terminus, in the river-like features themselves, and further inland regions (Fig. 8b). Any other assumed uniform fraction of overburden will similarly not account for spatial variations in pressure distribution.
(a) Difference in hydraulic head between water pressure assumed equal to 100% overburden pressure (p w/p i = 1.0) and that calculated in our main winter SHAKTI simulation. (b) Difference in hydraulic head between water pressure assumed equal to 80% overburden pressure (p w/p i = 0.8) and our main simulation. (c) Streamlines based on hydraulic potential flow routing with assumed 100% overburden pressure.
In the presence of steep topographic gradients, as encountered beneath Helheim Glacier, the deeply incised bed topography (and its reflection in surface topography) largely determines the locations of the main drainage paths. Recall that flow is driven by the hydraulic head gradient, and hydraulic head is composed of two components, pressure head and elevation head: h = p w/(ρ wg) + z b. In the case of localized canyons with bed elevation well below sea level, the elevation head in these troughs is sufficiently low to attract water from the surrounding areas, even when the water pressure is higher at the bottom of the canyons. As shown in Figure 8c, similar primary drainage pathways through the deep bed canyons emerge when using a simple routing calculation assuming water pressure equal to 100% overburden everywhere (compare to Re distribution in Fig. 2d). This is true using other uniform fractions of overburden as well, but the configuration of tributary drainage feeding the deep canyons differs, along with the resulting pressure distribution as shown in Figures 8a,b.
For purposes of identifying primary drainage paths, hydropotential routing provides a good estimate, assuming accurate bed topography, particularly in places with high-relief mountainous features. To investigate processes involving broad areas of the bed and drainage tributaries, however, spatially heterogeneous pressure distribution must be considered.
4.3 Winter priming of the drainage system
Our winter base-state main simulation highlights the fact that there is widespread subglacial drainage present year-round under Helheim Glacier, with at least a minimal system supported by basal meltwater generated by geothermal flux, frictional heat from sliding, and dissipated heat from the water flow, as considered here. At Sermeq Kujalleq (Store Glacier), a tidewater glacier in west Greenland, Cook and others (Reference Cook, Christoffersen, Todd, Slater and Chauché2020) also simulated an active winter subglacial drainage system, bolstered by winter discharge observations (Chauché and others, Reference Chauché2014). Our model results suggest that the drainage system at Helheim does not begin from a totally shut-down state at the initiation of melt each year, but retains some form through the winter. The ‘deflated’ drainage structure without surface meltwater input includes preferential pathways and large areas of low-transmissivity bed, primed to spring into more efficient drainage action with the delivery of surface-generated meltwater to the bed. With our winter simulation results in mind, the seasonal evolution of drainage efficiency and structure may not depend only on the spatio-temporal distribution of moulins and crevasses for meltwater inputs from the surface to the bed (or from englacial drainage of firn aquifers or other storage voids), but is likely also a function of the persistent winter base-state drainage structure. Therefore, we recommend considering an existing winter drainage system when interpreting observational data.
Previous work has explored the idea of winter priming beneath the Greenland Ice Sheet. In west Greenland, Chu and others (Reference Chu2016) found water stored on bed ridges in winter, which then flows to depressions and troughs in the melt season. Poinar and others (Reference Poinar, Dow and Andrews2019) simulated subglacial hydrology of an idealized Helheim-like glacier (without realistic bed topography) with year-round drainage of firn aquifer water to the bed, finding that increased water at the bed during winter facilitated more rapid and pronounced development of efficient channel networks in the melt season. The winter base state documented here would play a key role in that priming action, and shows that wintertime firn aquifer drainage may not be necessary in order to have year-round channelized structure at Helheim. This is largely because of the deeply incised canyons in the bed that preferentially pull water into river-like features even in the absence of meltwater inputs from the surface or englacial system, and this active winter channelized structure would likely be further enhanced by delayed meltwater drainage from the firn aquifer (Poinar and others, Reference Poinar2017). Recent work modeling Helheim subglacial hydrology by Stevens and others (Reference Stevens2022) also simulated active drainage in the near-terminus region of Helheim under winter conditions.
4.4 River-like winter features
SHAKTI employs a continuum description of subglacial geometry, without distinguishing between different drainage-system components. The way SHAKTI represents the geometry is through the subglacial gap height b, which is an average of the gap height over an entire element (i.e. generalizing earlier work on spatially lumped models such as Schoof, Reference Schoof2010; Brinkerhoff and others, Reference Brinkerhoff, Meyer, Bueler, Truffer and Bartholomaus2016). This means that SHAKTI does not resolve individual drainage channel geometry by calculating semi-circular cross-sectional area as in some other models (Hewitt, Reference Hewitt2013; Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013; Meyer and others, Reference Meyer, Fernandes, Creyts and Rice2016, Reference Meyer, Hutchinson and Rice2017; Felden and others, Reference Felden, Martin and Ng2023). The primary variable sought from subglacial hydrology models for ice dynamics calculations is the subglacial water pressure (or equivalently, effective pressure), which influences sliding velocity. The subglacial water pressure field is relatively smooth compared to small-scale geometric variations (e.g. in gap height at the bed) within a subglacial system. Without distinguishing different drainage modes with different evolution equations in each, SHAKTI is able to represent both distributed and channelized sub-systems naturally. With realistic bed topography incorporated, the winter simulation results presented above suggest the promise of SHAKTI in representing weakly connected sub-systems as well. The winter features reminiscent of broad channels have higher water pressure than their surroundings in this winter base state, but lower head, which is what drives the flux of water from the surrounding areas into these pathways due to the deeply incised bed. These river-like features will likely transition with seasonal meltwater input into even more efficient drainage channels. Stevens and others (Reference Stevens2022) also simulated substantial subglacial flow at Helheim with high-pressure channels persisting within 10 km of the terminus, even in winter, although that study considered a model domain focused on the near-terminus region without extending as far inland but employing a prescribed incoming flux from their upstream boundaries.
Note that the primary drainage pathways exist close to overburden pressure in winter. Recent work by Dow and others (Reference Dow, Ross, Jeofry, Siu and Siegert2022) simulated large, high-pressure channels extending far inland in Antarctica at close to overburden pressure, although their simulation results were sensitive to prescribed channel and sheet conductivities required by the GlaDS model. The existence of large, high-pressure channelized structures under thick ice is a relatively new idea, diverging from than the classical theory of low-pressure R-channels (Röthlisberger, Reference Röthlisberger1972). In some ways, interior regions of Greenland in the winter with no surface melt are similar to Antarctica. If Antarctica begins to experience more surface melt farther inland, this raises the possibility that it could potentially transition to become similar to Greenland in terms of seasonal hydrology behavior.
4.5 Role of basal stress in calculating frictional heat
Basal melt from frictional heat generated by sliding of the ice over its bed is potentially a dominant source of basal melt, as shown in Figure 4 and according to basal melt rates for Greenland calculated by Karlsson and others (Reference Karlsson2021), especially in fast-flowing glaciers like Helheim (which moves rapidly even in winter, with winter velocities exceeding 8,000 m yr−1) (Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017). The total basal melt rate we find for Helheim in the main simulation (6.3 m3 s−1, Table 4) is similar to the melt rate of 5.96 m3 s−1 calculated for Store Glacier in west Greenland by Cook and others (Reference Cook, Christoffersen, Todd, Slater and Chauché2020). This estimate was produced with a different model formulation and for a smaller, slower-flowing glacier than Helheim, but provides some corroboration for the theory that using driving stress or yield stress as proxies for basal stress may over-predict basal melt rates.
High localized basal melt rates on the order of 20 m yr−1 are calculated in our drivingstress and yieldstress simulations (Fig. S4), two orders of magnitude higher than the maximum melt rate in the nofrictionheat simulation without frictional heat from sliding, 1.0 m yr−1. In our original simulation with τ b = C 2Nu b, the maximum local melt rate is 14.4 m yr−1 (Fig. 2f). The plausibility of such high wintertime local melt rates is unresolved, as this rate of basal melt is inconsistent with most observations to date. Young and others (Reference Young2022) inferred basal melt rates of this order of magnitude in west Greenland, but in the context of a summer rain event, not a winter background melt rate. Greenland basal melt rates as calculated by Karlsson and others (Reference Karlsson2021) are typically < 0.25 m yr−1, which agrees better with our nofrictionheat simulation. In previous work on a Helheim-like idealized glacier, Poinar and others (Reference Poinar, Dow and Andrews2019) prescribed a uniform basal melt rate of 0.02 m yr−1 (based on thermo-mechanical modeling by Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012). In modeling subglacial hydrology at Helheim, Stevens and others (Reference Stevens2022) also prescribed a uniform basal melt rate of 0.0262 m yr−1 to represent geothermal and frictional heat contributions. The high melt rates in our simulations are extremely local; the average basal melt rate over a larger area is lower (mean of 0.097 m yr−1 over the entire model domain in our main simulation). With high localized basal melt rates included in an ice-dynamics model like ISSM, unusual features may result or the ice would need to compensate by flowing in to fill these melting ‘sinks’. This is an interesting topic to pursue in coupled simulations of hydrology and thermo-mechanical ice dynamics.
Given the disparity between the nofrictionheat simulation and the simulations using various formulations for basal shear stress, we demonstrate that frictional heat is an influential control on determining subglacial drainage regimes. This is not entirely satisfying from a modeling perspective. To add to the complexity, friction at the ice-bed interface may not be as straightforward as is frequently assumed, and heat may be generated deeper in the basal sediment as discussed by Hansen and Zoet (Reference Hansen and Zoet2022), which we do not include in our simulations. We must carefully consider what is likely to be a realistic drainage configuration at Helheim Glacier in the winter. Should we expect widespread high water pressure and clearly defined river-like pathways strongly influenced by the bed topography as we find in our main simulation? Or could there be sufficient heat generated by rapid sliding over the bed so that the water flow is actually more distributed and widespread through the main branches of the glacier without as distinct river-like features, as seen by using driving stress or yield stress for basal shear stress? Including frictional heat from sliding with these latter two methods leads to higher effective pressure (lower water pressure) over most of the interior domain (Figs. 6a and S5), contrary to our modeling target of reproducing high water pressure in winter. We are encouraged by the success of our main simulation that invokes drag coefficient and friction based on ISSM inversion in producing widespread winter water pressures and poorly connected regions of the bed.
Recent work by Stevens and others (Reference Stevens2022) to simulate subglacial hydrology at Helheim used a spatially uniform and temporally constant formulation to represent frictional heat over their domain, based on a value of basal shear stress of 60 kPa and a sliding velocity of 500 m yr−1. This approach does not capture the spatial variation in the contribution from frictional heat to basal melt, considering that much of the near-terminus region of Helheim has a persistent velocity of several km per year. We calculate basal stress lower than 60 kPa near the terminus and an order of magnitude higher (105 Pa) along the walls of the deep bed canyons (Fig. 5b).
Interestingly, the amount of meltwater produced in the domain is significantly greater when using driving stress as basal stress (drivingstress simulation) compared to the main simulation (Table 4). Over large scales, the basal shear stress should average to the driving stress, so these results imply that basal shear stress is concentrated in slower-flowing regions and along edges of fast-sliding areas (as seen in τ b from main, Fig. 5b). Employing driving stress as a proxy for basal stress in fast-flowing regions like the main branches of Helheim leads to an over-prediction of frictional heat. Admittedly, in simulations coupling effective pressure and ice velocity, the basal stress should match the driving stress well over large areas as the velocity adjusts to evolving effective pressure. In our stand-alone SHAKTI simulations assuming steady velocity, basal drag in the interior drifts from driving stress. This serves as excellent motivation for further work coupling hydrology and velocity. Future work employing a fully coupled model will provide improved insights on integrated mass-momentum-energy balances at the ice-water interface.
4.6 Seasonal meltwater inputs
To illustrate the capacity of the reduced SHAKTI model to handle transient meltwater inputs and the recovery period to winter conditions for interested readers, we also include results from a sample transient simulation with a duration of one year in the supplementary material (Figs. S6 and S7), applying distributed meltwater input (term i e→b in Equation (2)). The magnitude of meltwater applied to the bed is based on 2020 MERRA-2 data for a location in the main trunk of Helheim (GMAO, 2015), with hourly meltwater inputs smoothed over a period of 14 days. Beginning from the winter base state obtained through our main simulation, we supply seasonal meltwater inputs to the system at low elevations (surface elevation ≤900 m) to simulate the full annual 2020 cycle. Given that lower Helheim contains substantial crevassed regions, surface meltwater input does not necessarily reach the bed through point inputs such as moulins, as in western Greenland, and we have elected to approximate these low-elevation inputs as distributed evenly over the bed to represent widespread crevassing.
To maintain the focus of this work on winter base-state hydrology and the reduced SHAKTI model formulation, this simple transient simulation does not consider delayed drainage from the firn aquifer at higher elevations as in Poinar and others (Reference Poinar, Dow and Andrews2019) or lake drainages as in Stevens and others (Reference Stevens2022), but is included purely as an uncomplicated demonstration of the robustness of the reduced SHAKTI model to handle transient meltwater inputs and produce results within realistic ranges.
5. Conclusions
In this paper, we describe a reduced form of the SHAKTI subglacial hydrology model, retaining only the essential dynamics and parameterizations that do not involve poorly constrained parameters. We demonstrate the utility of SHAKTI by simulating the winter base-state drainage of Helheim Glacier in east Greenland. Like all models, SHAKTI is an approximation to a complex natural system, paving the way for large-scale coupled simulations. Stand-alone hydrology models such as SHAKTI are useful for understanding evolution of subglacial pressure and flow, but must be coupled to ice dynamics models to fully capture the reciprocal influences of effective pressure and glacier sliding. This coupling is an active area of ongoing work, and is beyond the scope of the present paper.
With the reduced SHAKTI model, we are able to: (a) reproduce widespread areas of high water pressure in winter using a continuum model, which are widely documented in field measurements and have been difficult to reproduce with subglacial hydrology models, and (b) demonstrate that hydraulic transmissivity, as calculated with the transitional flux formulation of SHAKTI, varies over several orders of magnitude within the domain, naturally representing poorly connected regions of the bed with a continuum approach.
Water pressure as a fraction of overburden varies substantially across the domain. While the locations of the main drainage pathways at Helheim are driven by the deeply incised topographic features in the bed and their corresponding effect in the ice surface slope, spatial pressure variations influence the overall flow configuration.
Frictional heat from sliding is the dominant source of basal melt rate over much of the domain, particularly in the interior, yielding high melt rates especially along the steep walls of bed incisions. The method used to calculate basal stress has a large influence on overall melt rate and the resulting drainage configuration and basal water flux. We find that using driving stress or yield stress as a proxy for basal stress both over-predict basal melt in fast-flowing areas compared to an approach that involves drag coefficients derived from inverse ice dynamics modeling. In summary, frictional heat from sliding is influential in the context of fast-flowing glaciers, and should be carefully considered in modeling subglacial hydrology and ice dynamics.
Finally, an important question for further research is to clarify the relationship between bed topography, basal melt rates, and maintenance of the subglacial hydrologic system.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2023.39
Acknowledgements
This study is part of a collaborative project funded by the Heising-Simons Foundation (grants # 2020-1911 (CRM, AS), # 2020-1910 (KP, JM), # 2020-1909 (WC), and #2021-3059 (MM)). SHAKTI is freely available for download as part of the Ice-sheet and Sea-level System Model at https://issm.jpl.nasa.gov/. Model output data and input scripts are available in a permanent Zenodo archive at https://doi.org/10.5281/zenodo.7019804. We are grateful for the Scientific Colour Maps developed by Crameri (Reference Crameri2021) and the Arctic Mapping Tools developed by Greene and others (Reference Greene, Gwyther and Blankenship2017). We thank Editor Ian Hewitt, along with reviewers Laura A. Stevens and Samuel Cook, for their helpful comments.
Open access
