Subglacial hydrology model overview

The subglacial hydrology model is based on the model described by Hoffman and Price13 for coupled distributed and channelized drainage. The model takes a continuum approach to subglacial drainage that captures the bulk behaviour of the subglacial drainage system at spatial scales relevant to ice dynamics (∼102 m) rather than attempting to model every basal conduit explicitly. The three modes of drainage, distributed, channelized and weakly connected, are modelled as separate components, each with appropriate physics, that are coupled together by exchanges of water driven by gradients in hydropotential between the systems. Distributed drainage is modelled as a continuous macroporous sheet in two dimensions on the primary model grid, while channelized drainage is represented by a single one-dimensional channel located along a line of grid cell edges of the primary model grid taking up no area within the primary model grid (Fig. 2b). Summaries of the previously described13, distributed and channelized drainage components are included in the Supplementary Methods for completeness. The weakly connected cavity system is represented as a subgrid element within the distributed system (Fig. 2b, inset). Within each grid cell of the model, a fraction of the bed, f w , is assumed to be covered by the weakly connected cavity system, with the remaining fraction, 1−f w , covered by the through-flowing distributed drainage system. Surface melt draining from the surface is delivered to the channelized system, which can exchange water with the distributed system. The distributed system can exchange water with both the channelized system and the weakly connected system, while the weakly connected system and the channelized system have no direct exchange. In all three systems, conduit space is assumed to be entirely filled with water at all times, a common approach in subglacial hydrology models (c.f. Schoof et al.51). All model parameters are listed in Supplementary Table 1 and are spatially uniform unless otherwise mentioned.

Weakly connected drainage model

The weakly connected component is implemented as a discontinuous reservoir that can exchange water with the surrounding distributed drainage system within each grid cell (Fig. 2b, inset). The weakly connected areas are prescribed to cover two-thirds of each grid cell, with the remainder covered by the distributed system. This fraction is chosen to capture the fact that all boreholes observed at drill site FOXX4,20 exhibited characteristics of hydraulic isolation from the active drainage system over most of the observing period while avoiding being overly prescriptive with the new component.

Using a macroporous sheet continuum formulation for the weakly connected system analogous to that used for the distributed system (Supplementary Methods), mass conservation of the water thickness in the weakly connected component (h w ) is described by the balance between locally generated melt (m w ) and exchange of water with the sheet (γ w ),

where ρ w is the density of water, A w is the area of the weakly connected system within each grid cell, defined as f w ΔxΔy.

Evolution of cavity space within the weakly connected component is the same as for the distributed system, a balance of cavity opening by sliding of the ice over bedrock bumps and close by creep of their ice roof:

where u b is the sliding velocity, A is the temperature-dependent rate factor for ice deformation, N w is the effective pressure in the weakly connected system, and h r and l r are parameters describing the height and wavelength, respectively, of bumps on the bed. The value for A for the basal layer (Supplementary Table 1) is approximated from results of borehole temperature and deformation observations and ice flow modelling performed by Ryser et al.14 where the basal layer is temperate and an enhancement factor of 2.5–4 (taken here as 2.7) is appropriate for ice from the late Wisconsin found at this depth.

Energy for local melting within the weakly connected component, m w , also matches that for the distributed system:

where G is the geothermal heat flux, is the basal traction vector and L is the latent heat of fusion of water.

For all components, hydraulic potential is defined as

where the subscript {d,c,w} indicates one of the distributed, channel, or weakly connected systems, respectively. g is the acceleration due to gravity (m2 s−1), z b is the bed elevation (m), p is water pressure (Pa) and p ice is the ice overburden pressure (Pa). For the distributed and channel systems the ice overburden pressure is calculated using a hydrostatic assumption, p ice =ρ i gH, where ρ i is the ice density (kg m−3) and H is the ice thickness (m). To allow the weakly connected system to attain water pressure greater than floatation, as observed, we include a simple representation for normal stress transfer34,36 in the weakly connected system by increasing the ice overburden pressure according to

We set r w to a constant value of 0.13, because this corresponds to the maximum value of borehole water pressure observed at our study site20, and thus suggests a limiting value when N w =0. Note that this simple representation does not include temporal variations in normal stress transfer.

Similar to the coupling between the channel and distributed system (Supplementary Methods), the weakly connected component is coupled to the distributed system by calculating a flux, γ w , between the surrounding distributed system and the weakly connected system within each grid cell based on a Darcy flow law:

where and are the hydraulic potential in the weakly connected and distributed systems, respectively, k 0w is the permeability between the weakly connected and distributed systems, P w represents the perimeter between the two systems within each grid cell and s is a characteristic spacing between the two systems.

While in reality the boundary between the two systems is likely to be complex, for simplicity in our parameterization we assume a simple geometry for the weakly connected patches within each grid cell for the purposes of generating reasonable values for P w and Δs. Specifically, we assume each weakly connected patch is a circle of radius 10 m (based on observations of differing connectivity on that scale in GrIS and elsewhere4,34,35,52) and set Δs=10 m. For the chosen values of the fractional area of the weakly connected system (f w =0.67) and grid spacing (Δx=Δy=200 m) we can calculate P w ≈5,355 m from basic geometry. While other choices could be made here, these details are not important from a practical standpoint without detailed knowledge of the bed, as these parameters, along with k 0w , are free variables in equation (6). Essentially, equation (6) parameterizes the exchange of water between the two systems primarily as a function of the difference in hydraulic potential and a permeability constant. As with the exchange between the distributed and channelized components (γ c ; Supplementary Methods), the exchange between the weakly connected and distributed systems can occur in either direction, with water always moving from the component with higher hydraulic potential to that with lower.

The permeability between the weakly connected and distributed systems, k 0w , is many times smaller than the permeability within the distributed system itself, k 0 (by ∼9 orders of magnitude in our model configuration; Supplementary Table 1), defining the weakly connected nature of this new component. To parameterize increases in connectivity between isolated regions of the bed and the active drainage system that are inferred to occur during summer, we allow k 0w to increase over the course of the summer by

where is a base value during winter, and k rate is a rate at which the permeability increases beginning at the start of summer, t s . We choose the value of k rate to reproduce observations of borehole water pressure (Fig. 4b). With our parameter choices, k 0w increases by a factor of ∼70 × over summer; the end of summer value of k 0w remains ∼8 orders of magnitude smaller than the permeability constant for the distributed system. While it changes in time, the permeability is spatially uniform. This is a simple parameterization for the increased connectivity in some boreholes observed during summer at site FOXX4,20 and at mountain glaciers34,35. It should be highlighted that there currently is little observational or physical basis by which to construct a governing relation for how permeability may evolve during periods of meltwater input, and improving on the simple linear relationship used here is a critical area for further research. Having explored a number of possible, more complicated ad hoc relations, we found that any relation that caused k 0w to increase substantially during summer generated similar qualitative behaviour.

Acknowledging that the weakly connected system may be a subset of distributed drainage but with very low permeability, an alternative implementation would be to directly model a distributed system with spatially varying permeability, avoiding the need for a new component. From a practical perspective, the subgrid parameterization of a third component used here is advantageous for two reasons. First, because GPS and satellite measurements of ice surface velocity indicate smooth velocity fields, it can be inferred that spatial variability in effective pressure and associated drainage conditions at the bed is at the length scale of an ice thickness or less. Representing such spatial heterogeneity at the grid scale would require a very fine grid, while the subgrid parameterization allows the weakly connected system to be represented at a coarser resolution, making this approach transferrable to large-scale ice sheet models. Second, explicitly modelling variable permeability of the distributed system would require additional assumptions and parameters about the spatial distribution of these variations.

Our parsimonious parameterization of the weakly connected system leaves room for additional complexity as empirical knowledge of the system increases. We have assumed in our model that all parts of the weakly connected system have changing permeability during summer (equation (7)), while borehole observations suggest that only some weakly connected cavities become ‘leakier’ during summer4. Similarly, some studies have observed weakly connected cavities becoming fully ‘connected’ during periods of high water pressure in the surrounding drainage system34,35, which would correspond to temporal changes in our f w parameter that we have kept steady in time (see Methods: model sensitivity to weakly connected area fraction for additional discussion). Additionally, the bedrock geometry for the weakly connected system may have different characteristics (h r , l r ) than for the distributed system. Though we acknowledge the inherent complexity of the subglacial system, we have chosen the simplest formulation that includes the dominant processes inferred to occur. Our parameterization of the weakly connected system is meant to represent the mean conditions of the bed and thus will not necessarily directly reflect unique measurements from specific boreholes.

Model setup

We generate a simplified model domain that is consistent with the basic geometry of our study site (Fig. 2). The domain represents a 100 km long sector of the GrIS from margin to lower accumulation zone, with our study site located 25 km inland from the ice margin. The domain is 5 km wide with periodic lateral boundaries to approximate the typical width of a supra- and sub-glacial catchment in this region of GrIS (estimated moulin density in this region is 0.2–0.25 km−1 (refs 1, 53)). The ice sheet geometry is a flat bed with a ‘plastic’ glacier shape51,54 assuming a constant basal shear stress of 105 Pa. This geometry provides an idealized but consistent setting with our study site (Fig. 2 and Supplementary Table 3). The subglacial hydrology model is applied for the entire domain, but model results are primarily analysed at the study site location (Fig. 2 and Supplementary Table 3); the larger domain is used only to generate realistic far-field constraints on the model solution at the study site. There is a zero subglacial water flux boundary condition for the distributed system at x=−100 km, and water pressure in the distributed system is assumed to be at atmospheric pressure at the ice sheet margin, x=0 km.

We choose this simplified geometry over the more complex geometry of the field site to simplify our interpretation of model results, given uncertainty in the true ice sheet geometry. Previous work in our study area has shown that complex topography can affect sliding and deformation over short distances14 and that ‘active’ and ‘passive’ subglacial regions might be affected by bed topography20. However, the two-dimensional basal topography of our study region is only known approximately. Only a handful of flight lines of airborne ice thickness measurements exist in close proximity to our study area, and for those that do, differences in ice thickness at flightline crossover points exceed 100 m in some locations. Additionally, the commonly used BedMachine product55, which uses mass conservation constraints to improve estimates of ice thickness between radar flight lines, does not cover this area. Previous subglacial hydrologic modelling has shown that valleys in the bedrock topography can concentrate subglacial drainage33, so the true bedrock topography is likely to play an important role in defining the large-scale drainage network and the initial rate of subglacial channelization. However, the location of moulins that input water to the basal drainage system is the key control on channel initiation and stability26,28,33. Because both our observations and simulations are focused on a site at a moulin, we expect the simplified model geometry to have little effect on the local subglacial drainage conditions at the study site.

To generate a model initial condition for the summer simulations, we spin up the coupled distributed and weakly connected systems to steady state, assuming no channelized drainage, to represent late winter conditions. We use the winter ice sliding speed to force the model and apply no meltwater forcing. Model parameters h r , l r and k 0 (Supplementary Table 1) are chosen to yield a hydraulic head in the distributed system at our study site of about 200 m below local floatation elevation, and is chosen to yield a hydraulic head in the weakly connected system of about 40 m above local floatation elevation.

Model forcing data

The subglacial hydrology model uses time series of two forcing fields, both of which are based on observations in the vicinity of the FOXX drill site. Surface melt input to the subglacial drainage system (ω in Supplementary Equation (5)) represents the volume flux of water drained to the bed through moulins and crevasses, and is the primary source of new water to the subglacial drainage system during summer (basal melting being of much smaller magnitude). Ice sliding speed (u b in equation (2) and Supplementary Equation (2)) controls the rate at which new cavity space opens in the distributed and weakly connected systems, and, less importantly, the frictional melt rate (equation (3) and Supplementary Equation (4)).

We base the melt forcing (Fig. 4a) off of 6-h average ablation rates4 measured using a pressure transducer installed below the surface connected to a surface reservoir56. Because the measured ablation rates do not have the precision necessary to directly generate a time-series with the required time step of 2-h, and the 6-h average rates reduce the amplitude of the diurnal cycle, we generate an idealized diurnal cycle by fitting a sine curve to the 6-h time-series where we adjust the amplitude on each day to maintain the measured 6-h average. The ablation sensor began to fail after the large melt event on days 229–230, so we estimate the ablation rate for days 231–240 using a linear fit between mean daily air temperature and ablation rate for the 10 days prior to the large melt event (Pearson correlation coefficient r=0.702, P<0.01). We scale the point ablation rate by the width of our model domain and apply it as a linear source term along the length of the channel model (ω). As discrete moulins can result in the generation of kinematic waves within our channel model, this linear distribution of melt input improves model stability, while reproducing synchronous changes in moulin water level as observed by Andrews et al.4 We apply a linear lapse rate such that runoff is zero at the equilibrium line altitude of 1,100 m. This provides plausible runoff rates for the entire domain, recognizing that model results will be most sensitive to runoff at the location of our study site which are well constrained. Though supraglacial storage on the GrIS is known to change over the melt season57, our melt forcing ignores supraglacial and englacial routing and storage processes that affect the timing and magnitude of melt delivery to the bed58,59, and is meant to be an idealized representation of the diurnal and seasonal variations in melt forcing at our study site.

The ice sliding speed forcing (Fig. 4a) comes from the results of Ryser et al.14 for site FOXX, which subtracts borehole-derived measurements of ice internal deformation from Global Positioning System-derived measurements of ice surface velocity to calculate basal slip. Four gaps of about a day are filled by averaging the diurnal cycles on either side of the gap and smoothing the edges to avoid any sharp transitions. Basal slip is 73% of motion during winter and, though speeds vary dramatically over summer, the contribution of deformation is roughly constant14. We apply this sliding forcing uniformly over our entire model domain. Though ice speed is not spatially uniform in reality, this simplifying assumption will not directly affect our conclusions because we only analyse the model results at our study site. The magnitude of basal traction, is held constant at 105 Pa for the entire simulation; though it should vary in reality, we do not have the information to provide a more accurate value, and, in any case, the results are not sensitive to this choice as it only affects frictional melting, which is orders of magnitude smaller than summer surface melt input.

Summer subglacial hydrology simulations

Using this spun-up state as an initial condition, we model the 2012 summer for days 150 (onset of summer speedup in the GPS record) to 250 (end of GPS and melt forcing observation time-series). For the summer simulation we add a single active channel along the centreline of the domain. While this precludes the formation of a network of channels (c.f. refs 26, 33), we expect this approach to resolve the dominant channel effects near our study site, as moulin location strongly controls channel nucleation28,33. The channel initial condition is an area of 10.0 m2 at the margin decreasing linearly to zero at 55 km inland, which yields an area of 5.45 m2 at the study site. This value was chosen to allow water pressure to drop below floatation and diurnal variations to develop within days of melt onset, as seen in the ice velocity record. The channel roughness parameter F is chosen to yield diurnal hydraulic head minima in late summer comparable to that measured in moulins4 and is broadly consistent with previous modelling efforts.

Hydraulic head observations

We compare model results to hydraulic head measured at site FOXX4,20 in 2012, where hydraulic head, d, relates to hydropotential as

Modelled hydraulic head in the channelized system is compared with hydraulic head measurements in moulin 3 and modelled hydraulic head in the weakly connected system to hydraulic head measurements in boreholes 4, 6 and 7. Because of the simplified model geometry and our desire to make the comparison of our model data general to both FOXX and neighbouring GULL study sites (which both exhibited similar borehole behaviour)4,20, there is no direct correspondence in ice thickness and bed elevation between the model and the observational study sites (Supplementary Table 3). Therefore, we plot both the model and measured hydraulic head using d f , the elevation corresponding to floatation pressure, as the vertical datum:

which is slightly different at each measurement site (see ref. 4, Extended Data Table 1) and in the model (Supplementary Table 3) due to modest differences in bed elevation and ice thickness.

Ice velocity calculations

Ice velocity for Fig. 3b,c is modelled using the thermomechanical, three-dimensional, first-order Stokes approximation60,61,62,63 momentum balance solver in the Community Ice Sheet Model v2.024 as described by Hoffman and Price13. The magnitude of basal traction, τ b , is defined by a physically based basal friction law for sliding over hard beds that allows for cavitation and bounded basal drag39,64,

where C is a Coulomb friction constant, λ max and m max are the wavelength (m) and maximum slope, respectively, of the dominant bedrock bumps, and n is the exponent in Glen’s flow law. is an addition to the original formulation added here to make basal traction calculated for our simplified model domain more realistic. Because our domain lacks the rough, heterogeneous topography of the real ice sheet where resistance to flow is likely to be concentrated in our study area20, we use to represent these missing ‘sticky’ spots40. We apply a constant value of across our study domain (Supplementary Table 1), chosen to reproduce a realistic range of model speeds across the range of effective pressure forcing applied.

With widespread cavitation (i.e., when effective pressure approaches 0 at high water pressure), the friction law becomes a Coulomb friction law of the form (moving to the left-hand side to clarify the form)

Alternatively, at large effective pressures (low water pressure) the friction law takes a power law form

The basal traction appropriate for ice dynamics is the integrated basal traction of both connected and isolated regions of the bed43,

where and are the basal traction in the distributed and weakly connected systems, respectively.

We approximate equation (13) by assuming the effective pressure used in equation (10) is an area-weighted average in each grid cell, N int , of the effective pressure in the distributed and weakly connected systems,

This simplification is justified by the fact that for the effective pressure and parameter values in the simulations, equation (10) remains primarily in the Coulomb friction law regime (equation (11)) where the use of equation (14) yields the appropriate description of basal traction, equation (13). (This approximation becomes increasingly inaccurate as the basal friction law transitions into the form of equation (12) where is not directly proportional to N).

Two standalone simulations of ice velocity are performed, both forced by effective pressure generated by the summer subglacial hydrology simulations. In the first simulation, equation (14) is calculated from N d and N w calculated in the summer subglacial hydrology simulation. This simulation demonstrates the impact of weakly connected drainage on ice dynamics (Fig. 3b). In the second simulation, the N w field is held steady while N d comes from the summer subglacial hydrology simulations (Fig. 3c). This is a control simulation that confirms the lack of seasonal hysteresis in the effective pressure–velocity relationship when the evolution of weakly connected drainage is not included in the model. Comparison of the velocity versus pressure relationship from these two simulations eliminates the possibility that the observed seasonal hysteresis is due to changing conditions in the distributed and/or channelized systems.

In both ice velocity simulations, the temperature-dependent rate factor, A, is a function of ice temperature40, and the vertical ice temperature profile in the model is taken equal to that measured at our study site14,22 using 11 uniform vertical levels. Based on ice deformation calculations from Ryser et al.14 we apply an enhancement factor to A of 2.7 × within the deepest modelled layer. Ice geometry is held steady for the duration of these summer simulations; ice velocity is simply calculated diagnostically at each time step based on the fixed geometry and changing basal boundary condition, which is a function of our modelled subglacial hydrological evolution. As for the subglacial hydrology model, the domain is periodic in the y-direction. Because we are only concerned with the ice speed at the study site, the model domain for the ice dynamics calculations is subset to span the area 10 km upstream and downstream of the study site to make the calculations less expensive. At the upstream and downstream boundaries we apply a vertically uniform Dirichlet boundary condition on the velocity field of 100 m a−1. We confirm that the velocity solution at the study site is independent of the choice of boundary condition value (as expected when the boundaries are far enough from the study site, >∼4–10 ice thicknesses65).

Parameters in equation (10) (Supplementary Table 1) are tuned to to yield model velocity of the observed magnitude, and the same values are used for both model versions. Specifically, and C are varied in combination to approximately match the observed diurnal range of surface speeds. The ratio is tuned to achieve the observed variation in the sensitivity of ice speed to effective pressure. This is assessed by an analysis of the slopes of the lines defining the minimum and maximum channel hydraulic head and ice surface speed on each day in Fig. 3. This slope represents the relationship between water pressure and sliding within a single day. A key feature in the observations (Fig. 3a) is a tendency for lines restricted to lower hydraulic head values to have flatter slopes than those at higher hydraulic head values (lines to the left tend to be more horizontal). This behaviour is the result of a reduction in sliding sensitivity to water pressure at low water pressures (high effective pressures). It is predicted by theory and is a key feature of the basal friction law used (equation (10)). At large effective pressure when subglacial cavities are small, sliding is controlled by regelation and enhanced creep38,39,40. This insensitivity of velocity to effective pressure at large effective pressures has been observed in mountain glaciers41 and can be clearly seen in high temporal resolution data at our study site (see ref. 4; Extended Data Fig. 4b). To ensure the models are calibrated to correctly represent this changing sensitivity of sliding to effective pressure, we calculate a linear regression between the minimum channel hydraulic head and the slope of the channel hydraulic head/ice speed relationship on each day. Restricting the regression to the range of hydraulic head minimum values in common between all three data sets (<125 m), we confirm that the observations and both model versions have a similar sensitivity to changes in effective pressure for the parameter values used (Supplementary Fig. 1). The differing slopes of the lines in Fig. 3c represent this variable sensitivity of sliding to effective pressure and are what cause the modest seasonal-scale changes exhibited by the model with a static weakly connected system.

Model sensitivity to weakly connected area fraction

While a complete analysis of the sensitivity of all parameters used in the models is beyond the scope of this study, we assess how the study’s main conclusions are affected by the choice of key parameter, f w , the area fraction of the weakly connected system, and the possibility that it could change in time.

In the main text, we present results for f w =0.67. Because a total of six boreholes at two different drill sites exhibited out-of-phase behaviour for the majority of both summers measured4,20, we assume the fraction must be large, but we want to avoid being overly prescriptive with the new component. In fact, an upper bound for f w can be found based on the observations4,20 that hydraulic head in the moulin system has typical diurnal variations of about 20% of overburden pressure, while diurnal variations in the boreholes are about 2%. Because the observed surface velocity is precisely in-phase with the moulin pressure variations4 that are driving the system, we can assume that the area-weighted diurnal variations in the distributed system connected to the moulin must be larger than the area-weighted variations in the weakly connected system:

This provides an upper bound of f w <0.91. Note that if diurnal variations in the moulin are in actuality larger than those across representative regions of the distributed system connected to it (as might be expected for a diffusive pressure wave32), then the upper bound for f w should be smaller than the calculated value.

As anticipated from equation (15), if we prescribe f w =0.90 within the model (and retuning the parameters in the basal friction law), diurnal variations in ice speed are almost entirely absent (Supplementary Fig. 5). This is because the area-weighted diurnal variations of the weakly connected system roughly cancel out the area-weighted diurnal variations in the distributed system because the two systems are out of phase from each other. The few days with substantial diurnal range in ice speed seen in these model results occur when conditions within the model temporarily shift the phasing of diurnal variations in the weakly connected system. Note that ice speed still drops over the summer when the weakly connected system is allowed to drain (Supplementary Fig. 5b), but overall the results differ markedly from the observations and are largely unphysical.

The lower bound on f w is less constrained than its upper bound, but based on the weakly connected behaviour of all boreholes in the study area, we assume f w =0.50 forms a reasonable lower constraint. Model results with that value (Supplementary Fig. 6) are somewhat similar to the baseline value of f w =0.67 used in the main text (Fig. 3), but the ice speed does not drop as substantially; with the weakly connected system covering a smaller fraction of the bed, changes there have less impact on the integrated basal traction. Noting that the modelled changes in the weakly connected system have been constrained by the borehole observations, we find that f w values of ∼0.60–0.80 can give results that provide a reasonable match to the ice speed measurements, allowing for modest adjustments to the other parameters in the model.

In addition to assessing the baseline value of f w , we also consider the possibility that f w could change during the summer. Certainly, rather than existing areas of weakly connected drainage becoming more strongly connected to the active drainage system, a reasonable hypothesis would be a change in the area fraction of weakly connected regions. We test this by an additional model run where evolution within the weakly connected system occurs by f w declining (Supplementary Fig. 2) or increasing (Supplementary Fig. 3) rather than changes to the permeability. These parameterizations do a worse job at reproducing the observed behaviour. Of these two runs, the situation where f w declines is the more observationally supported change—observations on mountain glaciers have shown that ‘isolated’ boreholes can become connected during periods of high water pressure in the active drainage system34,35. However, in the run where f w declines, little seasonal evolution in the ice speed occurs (Supplementary Fig. 2b). Instead, the diurnal range in ice speed increases as the summer progresses, an effect that is not seen in the observations. Similarly, the primary effect of increasing f w during summer (Supplementary Fig. 3b) is a decrease in the diurnal range of ice speed. While we cannot rule out the possibility that the area fraction of the weakly connected system changes during summer, our model results indicate it is not the primary mechanism causing ice speed to drop.

Data and code availability

Model code is development code based off of the Community Ice Sheet Model (CISM) version 2.0.4 (http://oceans11.lanl.gov/cism/). Model code, processing scripts, input datasets and model output are all available on request from the corresponding author.