Water pressure fluctuations control variability in sediment flux and slip dynamics beneath glaciers and ice streams

Rapid ice loss is facilitated by sliding over beds consisting of reworked sediments and erosional products, commonly referred to as till. The dynamic interplay between ice and till reshapes the bed, creating landforms preserved from past glaciations. Leveraging the imprint left by past glaciations as constraints for projecting future deglaciation is hindered by our incomplete understanding of evolving basal slip. Here, we develop a continuum model of water-saturated, cohesive till to quantify the interplay between meltwater percolation and till mobilization that governs changes in the depth of basal slip under fast-moving ice. Our model explains the puzzling variability of observed slip depths by relating localized till deformation to perturbations in pore-water pressure. It demonstrates that variable slip depth is an inherent property of the ice-meltwater-till system, which could help understand why some paleo-landforms like grounding-zone wedges appear to have formed quickly relative to current till-transport rates. The disparity in predicted and observed till transport rates beneath flowing ice is partially explained by water pressure fluctuations which alter the rate and depth of slip in the sediments, according to a coupled ice-meltwater-till continuum model

T he Intergovernmental Panel on Climate Change ascribes the largest uncertainty in future sea level rise to the ice dynamics of our two ice sheets, Greenland and Antarctica 1 . Most of this uncertainty arises from specifying meaningful boundary conditions at both the base of the land ice and its seaward margins 2 . While ice-ocean interactions have received considerable attention over recent years [3][4][5] , our ability to predict the evolution of basal conditions underneath the ice sheets remains limited.
The debate over the evolution of basal conditions is closely tied to the question of where basal slip is accommodated. The fast speed of the ice leaves little doubt that motion is facilitated primarily by localized basal slip rather than distributed deformation inside the ice, but observations of where slip occurs provide conflicting evidence. Field observations range from basal slip occurring almost exclusively at the ice-bed interface 6,7 , to partial slip at the ice-bed interface combined with shallow deformation within the upper portion of the till 8,9 , and deep basal slip meters below the ice-bed interface 10,11 . Interestingly, the depth of basal slip can vary by several meters even at the same field site as in the case of Black Rapids Glacier, Alaska 9,10 . Observations hence suggest that the depth of basal slip is a dynamic variable that evolves in time.
The mechanical coupling between the ice and bed determines if till is mobilized beneath the ice-bed interface [12][13][14][15][16] . Coupling could emerge, for example, from protruding clasts lodged in the basal ice that plough through the uppermost bed sediment [17][18][19][20][21][22] , and from regelation infiltration incorporating the upper till bed into the basal ice through freeze-on 23 . The depth of basal slip in the till bed determines both the rate and the volume of till transport, and directly influences the pace at which a glacier or ice stream reshapes its bed [24][25][26][27] . The landforms emerging from this reshaping process alter ice motion and till deformation by locally changing the bed slopes, rerouting subglacial water [28][29][30] and, potentially, creating barriers to ice transport 31 . Leveraging the information imprinted in till and landforms deposited from past glaciations 32,33 to inform models of current and future deglaciation, requires understanding the feedback between ice motion, till deformation and water storage and drainage. The seemingly inconsistent observations of basal deformation across space and time mean that our understanding of the dominant subglacial processes is incomplete, shedding doubt on the reliability of predictive models.
The goal of this paper is to advance our fundamental understanding of the physical processes that contribute to the dynamics of subglacial beds as reflected in existing, seemingly inconsistent observations of basal deformation. We hypothesize that the variability in observational evidence reflects different dynamic regimes arising from the complex interplay between ice, meltwater and till in the subglacial environment. We test this hypothesis by deriving a continuum model for water-saturated subglacial till that is consistent with laboratory measurements of a near-plastic till rheology 16,[34][35][36] . Our model builds on recent advances in continuum descriptions of dense, granular behavior by Henann and Kamrin 37 and integrates it with progress in understanding the coupling of water flux through granular media at the scale of individual grains [38][39][40] .
The granular mechanics of cohesive, water-satured till Over the past two decades, many granular mechanics models have shifted towards statistical approaches 41 and away from classical solid mechanics 42 . The main advance provided by this line of thinking is the ability to represent non-local effects, which represent deformation in places where the granular medium is not strictly at the yield stress but deforms due to interactions with failing patches nearby. These non-local interactions dominate deformation in granular shear zones [43][44][45] such as those under fast-moving ice. In the non-local granular fluidity (NGF) model by Henann and Kamrin 37 , a fluidity field variable accounts for the non-local effects on deformation. Microscopically, fluidity is related to velocity fluctuations in the granular medium 46 .
The original non-local fluidity model accurately describes the strain distribution in a variety of experimental settings 37 , but it assumes a dry, cohesionless granular medium. In the subglacial context, in-situ measurements document significant temporal variations in both water pressure and shear strain rates 6,8,10,47 . These variations result from both internal till dynamics 48 and external water input 49 , highlighting that meltwater is a crucial component of the overall dynamics. In the current study, we capture the interplay between till strength and water percolation by incorporating effective stress and pore-pressure diffusion that depends on porosity and permeability 38,39 . We also consider cohesion in subglacial tills, which tends to increase with clay content 50 , by adding a strength contribution from cohesion to the model. We refer to our model as the cohesive non-local granular fluidity model with pore fluid (CNGF-PF). Contrary to the intense computational costs of our previous particle-scale model (e.g. 51,52 ), the CNGF-PF model is a continuum model that is sufficiently efficient to allow coupling to numerical ice-dynamics models. Following 37 , we assume that till is in the critical state throughout the domain, implying that average porosity does not change as a function of granular deformation. The critical state assumption is likely justified, particularly underneath fast-moving ice 16 .
To test the ability of our model to reproduce the range of observed behavior in the field, we shear the modeled till with idealized boundary conditions that represent end members of the ice-bed coupling (Fig. 1a). Speed-controlled conditions are typical for laboratory shear experiments, and represent a subglacial setting where changes in bed friction do not immediately alter the ice flow velocity because the force balance includes other prominent terms such as topography or shear stress in the lateral margins. Specifically, the interfacial shear velocity at the base may vary over time, but is not dependent on the local basal resistance contributed by subglacial friction. An example of this limit is Storglaciären Glacier, Sweden, where subglacial strain rates and surface velocities are largely uncorrelated, at least at the place and time of measurement by Hooke et al. 8 . Similarly, the slip speed of several outlet glaciers in Greenland show little sensitivity to variations in meltwater production 53 . In contrast, stresscontrolled simulations approximate conditions where ice flow velocity responds to changes in subglacial strain rate directly. Whillans Ice Stream, West Antarctica is an example of this setting, where a low surface slope and low driving stress results in stick-slip movement 54 . Evidently, most real glacier settings shift along the spectrum spanned between the two limits, depending on how important basal friction is to the overall stress balance.

Results
Consistency with laboratory measurements and granular-scale simulations. Before applying our model to field data, we validate it against laboratory measurements of the behavior of glacial till and DEM models. The validation procedure relies on the effects of two non-dimensional parameters of the CNGF-PF model: A scales the spatial extent of non-locality, and b scales the weak rate dependence on friction beyond yield. Natural glaciers systems show extremely variable strain rates of up to~5 × 10 3 a −1 . When studied over five orders of strain-rate magnitude (Fig. 1b), some tested natural tills show slight rate weakening (Storglaciaren till, Two Rivers till, WIS till, Lower Cromer till), while others are slightly rate strengthening (Caesar till, Crowden till). Much like the natural tills, our simulated till is effectively rate-independent over most of the range. At extreme shear-strain rates (Fig. 1b), higher values of b provide larger frictional resistance, leading to a slight hardening. Using published, rather than fitted, values of internal friction μ s and cohesion C, our model reproduces the Mohr-Coulomb behavior of tills exactly [34][35][36] (Fig. 1c).
Laboratory constraints for the dimensionless non-locality parameter A are not available, because we are not aware of any experiments that analyze the effect of normal stress on the strain distribution in till. Instead, we compare the modeled strain distribution with discrete-element method results from 51 . By inserting relevant material parameters for grain size, friction, stress, and shear velocity (DEM parameters in Supplementary Table 1), our model approximates the strain distribution well when A = 0.5 ( Fig. 1d), consistent with models of deforming glass beads that used A = 0.48 37 . Both models show that sediment transport is pressure dependent, with low effective normal stresses producing shallow deformation, and high effective normal stresses deepening the deformation. Deep deformation leads to a higher sediment transport, because the entire sediment column above the shear zone is transported downstream along with the ice movement. However, the DEM results required more than two months of computational time on a modern GPU, whereas the continuum model is completed in a fraction of a second, albeit without details of individual particle kinematics and dynamical adjustment towards the critical state.
We emphasize that a local-only granular model, such as μ(I) 55,56 , does not include the effects of sediment texture (size, shape, and arrangement of grains) on shear-zone decay between actively flowing and stagnant parts of the bed. Comparing Fig. 1d (NGF) with Supplementary Fig. 1 (μ(I)) demonstrates that the non-locality significantly contributes to shear strain distribution in a set of example shear tests, highlighting the value of using a non-local rheology for the following analysis of subglacial deformation patterns.
Meltwater pulses can lead to deep slip in speed-controlled settings. To test the hypothesis that meltwater input can alter the depth at which basal slip is concentrated, we impose daily sinusoidal variations in water pressure at the top of the model domain and simulate pore-pressure diffusion and the resultant shear dynamics over seven days ( Fig. 2a and Supplementary Fig. 2). The effective normal stress generally increases with depth due to increasing weight of the overlying sediment. However, the timevariable input in meltwater at the ice-bed interface and its downward diffusion modulate the depth-dependence of deformation as well. Due to diffusion, water pressure perturbations decay exponentially with depth, and travel with a phase shift (Fig. 2d inset). As a consequence, the distribution of water pressure can create minima in effective normal stress at significant depth below the ice-bed interface, such that the effective stress profile locally reverses (Fig. 2d).
Under speed-controlled conditions (Fig. 2b, e), simulations show highly time-variable deformation patterns throughout the till. Granular failure and maximum shear deformation occur where the effective normal stress is at its minimum, resulting in a plug-like motion of the sediment layer above the slip zone during reversal of the effective stress at depth. Deep deformation and plug-like till transport occurs predominantly when the water pressure at the ice-bed interface is close to its minimum and rising.
Under the same water-pressure forcing, the stress-controlled conditions (Fig. 2c, f) produce a very different, dynamical pattern of deformation. Shear occurs almost exclusively in the uppermost part of the modeled till layer. The cyclic effective stress in combination with the frictional properties of the till lead to stickslip behavior, which is associated with dramatic variations in the shear speed (Fig. 2c). When water pressure at the top decreases and effective stress increases, shear-strain rates decrease from the top and are maximal just below the ice-bed interface. When the effective stress at the top further increases, the shear stress is Till flux is controlled by the coupling at the ice-till interface. In Fig. 3, we replot the simulations from Fig. 2 as time stacks to clarify the relationship between water pressure and stress, strain rate, and till flux. Under speed-controlled conditions, the shear stress varies in an approximately linear manner as predicted for Mohr-Coulomb materials 42 , but with hysteresis at high effective normal stresses (Fig. 3a). Shear stress is higher during deep shear because the lithostatic contribution increases effective normal stress and shear strength at depth. When driven by stresscontrolled conditions, the simulated water-sediment system shows stick-slip behavior with locally very high slip speeds and strong hysteresis (Fig. 3d). We emphasize that the specific values for the slip speed and till flux are not realistic for a glacial setting. Far-field stresses such as lateral support from the margins 57 , which are not included in our model, are expected to restrain the high till shear speed, and real ice-bed interfaces can transition from a stresscontrolled to a speed-controlled configuration during slip.
Under stress-controlled conditions, deformation occurs at shallow depth with rapid slip when ice-bed interface when water pressures are at their highest magnitude (Figs. 2c, f and 3c). The trends in when and how much till transport occurs are therefore drastically different for the two cases. The largest sediment transport under speed-controlled conditions occurs with transiently low water pressures and high effective normal-stresses at the ice-bed interface (Fig. 3b). Instead, the majority of till transport under stress-controlled conditions occurs as shallow deformation during rapid slip events and high water pressures (Fig. 3d).
Deep slip depends on water-pressure variations and till permeability. While the previous two sections shed light on the effect of meltwater pulses and the importance of the boundary condition imposed at the till-ice interface, the occurrence of deep slip also depends sensitively on the hydraulic permeability of the  till layer which can vary by orders of magnitude. Here, we derive an analytical solution that can provide a quick order-ofmagnitude estimate for the depth at which deep slip is expected to occur as a function of till permeability and amplitude and frequency of the water pressure variations.
The deepest deformation under speed-controlled conditions depend on the diffusion of pore-pressure perturbations away from the ice-bed interface. The length-scale of pressureperturbation decay due to diffusion is characterized by the skin depth d s [m]. Assuming that fluid and hydraulic skeleton properties are constant and the layer is sufficiently thick: where k [m 2 ] is the hydraulic permeability, α [Pa −1 ] is the grainskeleton compressibility, β f [Pa −1 ] is the compressibility of water, ϕ [-] is the porosity, η f [Pa s] is the dynamic viscosity of water, and f [s −1 ] is the frequency of the pore-pressure variations at the ice-bed interface. Figure 4a shows the skin depth for water at 0 ∘ C under a range of permeabilities and water-pressure forcing frequencies. Importantly, the deformation pattern depends sensitively on the distribution of effective stress and hence on the pressureperturbation amplitude A f [m], which means that the skin depth alone is insufficient to infer whether deep deformation occurs. We derive an analytical solution for diffusive pressure perturbation to find the largest depth beneath the ice-bed interface where the effective stress is minimal over the course of a pressureperturbation cycle (see Supplementary Note 2 for the full derivation). As the shear zone under speed-controlled conditions follows the effective normal stress minimum, this depth corresponds to the deepest expected position of the shear zone midpoint, z 0 , during a cycle of water-pressure forcing at the icebed interface: where ρ s and ρ f [kg m −3 ] is the till and water density, respectively, and G [m s −2 ] is the gravitational acceleration. The prediction of the analytical solution matches our simulations well (depth marked by dashed horizontal line in Fig. 2). We plot solutions to Eq. (2) for various combinations of the water-pressure perturbation amplitude (A f ), perturbation frequency (f), and permeability (k) in Fig. 4b and c. For a given pressure-perturbation amplitude, deep deformation does not occur above a threshold value for the skin depth. Beneath this threshold, deformation quickly deepens to a substantial depth before gradually shallowing at lower skin depth values. Our analysis hence suggests that deep deformation can occur in different tills. In high-permeability tills, deep deformation and till transport is facilitated by rapid changes in water pressure, while transport of low-permeability tills are sensitive to longer, seasonal pressure variations. Larger water-pressure magnitudes facilitate deep deformation at even higher till permeabilities and slower pressure perturbations.
Comparison to field observations Our model results demonstrate that deep slip in till occurs when remnant high water pressures at depth overcome the effective lithostatic stress gradient. The depth of slip for deforming beds is hence dependent on several evolving parameters including ice thickness, variability and magnitude in meltwater influx, and till properties. However, the key factor is the coupling at the ice-till interface or, in other words, the degree to which the slip speed responds to in-situ changes in the interfacial shear stress. Interestingly, in natural systems, the degree of resistance provided by the interfacial shear stress varies widely from being an almost negligible contribution to being the dominant term in the force balance 57 . The coupling could also change and, in fact, is likely to change over the course of a deglaciation or a glacial cycle.
The ice streams of the Siple Coast, West Antarctica, provide an interesting example in this context. While observations indicate that present till fluxes are small (<40 m 3 m −1 a −158 ), large sedimentary grounding-zone wedges on the continental shelf appear to have formed rapidly, suggesting a significantly greater till flux (1.4 × 10 3 m 3 m −1 a −1 ) during and immediately after the last glacial maximum 59 . We conceptualize the contemporary Whillans Ice Stream at the Siple Coast as a stress-controlled setting, which is consistent with its current stick-slip behavior 54 . Stresscontrolled subglacial deformation results in shallow shear below the ice-bed interface (Figs. 1d and 2c). The shallow till deformation at the ice-till interface produces relatively small till fluxes, which is consistent with modern observations 58 .
During the last glacial maximum, however, the ice thickness in the Ross Sea was larger than today and likely characterized by greater surface slopes. Both factors increase driving stress, suggesting that lateral stresses in the margin could have played a more important role in the force balance of the ice stream than they do today, at least given comparable till properties underneath. Therefore, the ice-till interface during the last glacial maximum could have resembled speed-controlled conditions more closely. A speed-controlled simulation produces shallow shear if the subglacial water pressure is steady (till parameters in Supplementary Table 1). Daily or weekly variations in ice-bed interface water pressures are too frequent to create significant perturbations in the low-permeability till underneath Whillans Ice Stream, but variations of approximately 50 kPa on the Fig. 4 Analytical solutions to skin depth and maximum deformation depth from a sinusoidal water-pressure forcing with amplitude A f at the ice-bed interface. a Skin depth of pore-pressure fluctuations (Eq. (1)) with forcing frequencies ranging from yearly to hourly periods. The permeability spans values common for tills. b Depth of maximum deformation with a water-pressure forcing of A f = 10 kPa. c Same as b but with A f = 100 kPa. monthly, annual, or decadal time scale can transiently deepen shear deformation to 0.2, 0.6, and 1.3 m, respectively, greatly enhancing till transport. Such longer-term variations, mainly driven from internal dynamics, have been observed both locally by in-situ borehole measurements 60 , and inferred over larger areas by remote sensing 61 , and are to be expected in particular during ice advance and retreat.
In our idealized analysis, we do not explicitly account for icebed interface physics of ploughing and regelation infiltration that can influence till flux during shallow slip [20][21][22] . However, we demonstrate that hydraulic variability has the potential to greatly increase subglacial till flux through deep seated deformation, which could contribute to rapid changes in sediment transport and the pace of landform development at the grounding zone. While we focus on depth-variability in deformation, horizontal variations are likely to arise as well 62 , particularly if subglacial water transport is heterogeneously distributed in space. The interplay between water pressure fluctuations and sediment mobilization together with a shift in the relative importance of basal strength could reconcile the conflicting observations of small till flux today with drastically larger fluxes under past glacial configurations.
Another interesting and observationally better constrained field site is Black Rapids Glacier, a contemporary mountain glacier in the central Alaska Range. A transition in the coupling at the icetill interface was observed in between field campaigns at Black Rapids Glacier that were only a few years apart. The till bed of Black Rapids Glacier was probed by in-situ borehole measurements repeatedly 9,10 . During the earlier field campaign in 1997 10 , basal hydrology and deformation did not significantly correlate with ice-surface displacement, suggesting that the subglacial bed can be approximated by speed-controlled conditions.
Assuming a diffusivity of 1.5 × 10 −5 m 2 s −1 , water-pressure perturbations of approximately 1 MPa, and a forcing period of 1 month 10 , we estimate the maximum deformation depth with Eq. (2) to be on the order of z 0 ¼ 7:2 m. This result is in broad agreement with observations that ice slip resulted from till deformation at an undetermined depth of more than 2 m below the ice-bed interface 10 . Our model suggests that significant volumes of till would be transported in this limit, because the entire till layer above the deep slip zone is mobilized. In contrast, later measurements at the same location on Black Rapids Glacier in 2002 9 showed that sliding occurred partially at the ice-till interface itself and partly at shallow depth within the till layer. The glacier geometry and flow remained largely unchanged between the two field campaigns, but the recorded water-pressure variations were smaller in 2002. Our model framework reconciles these seemingly contradictory observations, because it clarifies the role of varying pore-water pressure amplitude in shifting the depth of basal slip away from ice-till interface and deep into the till layer itself.
Our analysis highlights the value of continuous instrumentation of field sites 9,10,47,60 , because it suggests that the observed, large variability is an inherent property of the ice-meltwater-till system. As a consequence, measurements performed at a single point in time do not always provide sufficient constraints to assess the evolution of basal slip and till transport over longer time scales. Similarly, our results highlight the importance of not reducing the dynamics of the subglacial environment only to subglacial hydrology. While subglacial hydrology is important 30 , the effect it has on ice speed and the depth of basal slip depends sensitively on the coupling at the ice-till interface. This insight could be relevant for understanding the distinct but variable response of glaciers to seasonal water input 53 and to climate forcing more generally.
We emphasize that the motivation for providing a comparison between a highly idealized model and field observations is to demonstrate the explanatory potential of the model we have derived. There is no doubt that any actual field setting is characterized by significantly more complexity than considered here. Nonetheless, we argue that it is valuable to test how much of the observed variability can be reproduced within a relatively simple, unifying framework, where the complex coupling between ice, till and water is not confounded by additional factors such as topography, ice-ocean interactions, and hydrological processes to only name a few.
To conclude, we note that a key, general outcome of our analysis is the understanding and recognition that contrasting field observations of basal behavior do not necessarily indicate different basal physics. Instead, variability in till dynamics and flux could reflect the natural temporal and spatial variability in the effective stress conditions within the till due to changes in ice overburden and meltwater percolation. Importantly, our modeling framework provides an opportunity to further refine these insights, because it is computationally cheap and adjustable to various grid setups, making it conducive to integration with largescale ice sheet models.

Methods
Model formulation. Subglacial tills are saturated by meltwater and often have significant cohesion from clay minerals. To generalize the existing non-local, granular fluidity (NGF) model by 37 to subglacial till (original model detail listed in Supplementary Note 1), we allow pore pressure to perturb the stress distribution in the sediment, and add the contribution of cohesion to shear strength. The extended model is named the cohesive NGF pore fluid (CNGF-PF) model.
Typical shear-strain magnitudes in subglacial tills are >>1 16 , so we neglect the elastic contribution to shear strain rate (_ γ ¼ _ γ p [s −1 ]). The transient evolution of pore-fluid pressure (p f [Pa]) is modeled through Darcian pressure diffusion 38,39 : When comparing our model to published results of hydraulic diffusivity D [m 2 s −1 ], the right-hand side of Eq. (3) is D∇ 2 p f .
In the granular phase, we use the effective normal stress σ 0 ¼ σ n À p f [Pa] instead of the dry normal stress σ n . We further add a cohesion contribution C [Pa] to the shear strength in the cooperativity length, ξ, in 37 : and in the local fluidity term, g local : g local ðμ; σ 0 Þ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi d 2 σ 0 =ρ s q ððμ À C=σ 0 Þ À μ s Þ=ðbμÞ if μ À C=σ 0 > μ s ; and 0 i f μ À C=σ 0 ≤ μ s : The fluidity field is solved in the same manner as the original model 37 .
Simulation setup. We apply the model in a one-dimensional setup with simple shear (Fig. 1a) of length L z . The upper normal stress (i.e. the ice weight, σ n (z = L z )) is constant, and normal stress linearly increases with depth due to sediment weight. In the stress-controlled experiments, the upper boundary (i.e. the ice-bed interface) applies a fixed shear stress τ on the till. In speed-controlled experiments, the upper boundary moves with a fixed shear velocity v x (z = L z ). For all simulations, the lower boundary condition for the granular phase is no slip (v x (L = 0) = 0). Effective normal stress (σ 0 ¼ σ n À p f ) varies if the water pressure p f changes (Eq. (3)). For the water-pressure solver, the top pressure (p f (z = L z )) is either constant or varies through time. For the experiments with variable water pressure, we apply a water-pressure forcing amplitude of A f = 80 kPa that modulates effective normal stress at the top around 100 kPa. At the start of each simulation, water pressure is set to follow the hydrostatic gradient at the lower boundary (dp f /dz(z = 0) = ρ f G). Unless noted, we use the value 0.94 for the dimensionless rate-dependence parameter b, which is empirically constrained from laboratory experiments on glass beads 63 . We solve the model components iteratively as detailed in the Supplementary Note 1 with specific parameter values listed in Supplementary Table 1. A parameter sensitivity test is provided in Supplementary Note 3 and Supplementary  Fig. 3.

Data availability
No datasets were generated or analyzed during the current study.

Code availability
The grain-water model is written in C and is available under free-software licensing at https://src.adamsgaard.dk/cngf-pf (permanently archived at doi:10.5281/ zenodo.4106566). Run "man cngf-pf" after installing for usage information. The results and figures in this paper can be reproduced by following the instructions in the experiment repository for this publication, available at https://src.adamsgaard.dk/cngfpf-exp1 (permanently archived at doi:10.5281/zenodo.4106575).