A slip law for hard-bedded glaciers derived from observed bed topography

Numerical simulations of glacier slip over diverse hard beds indicate Coulomb slip behavior like that of deformable beds.


INTRODUCTION
Major contributions of ice discharge to the oceans that result largely from rapid slip of glaciers over their beds are accelerating mass losses from ice sheets and associated sea-level rise (1,2). Both slow-moving (3) and fast-moving (4-6) parts of ice sheets can rest wholly or in part on hard beds that behave rigidly so that ice slips over them rather than deforming them. Where fast-flowing ice rests mostly on deformable till ("soft" beds), hard-bedded zones can have disproportionately large effects on modeled ice losses (7). Moreover, glacier slip over hard beds causes erosion of bedrock that modulates uplift rates in mountain belts (8) and rates of chemical weathering that alter climate (9).
Models of glacier flow in response to climate forcing (2,10) and of orogen-scale glacial erosion and its effects (11,12) therefore require a slip law for hard beds that relates drag at the bed,  b , to slip velocity (13). Modeling studies indicate acute sensitivity of glacier response to the form of this relationship (2,(14)(15)(16), which is affected by water-filled cavities that persist in the lees of bumps on the bed where ice pressure on the bed is low (17). The cavity water pressure, subtracted from the ice-overburden pressure, yields an effective pressure, N, that also contributes to the balance of forces on the bed. Quasi-static equilibrium requires that for two-dimensional (2D) beds, the ratio,  b /N, cannot exceed the maximum slope, m max , of the upglacier (stoss) sides of bed bumps (18).
Process-based models for hard-bedded slip with cavities indicate slip laws with a common form: With increasing slip velocity or decreasing N,  b /N increases, peaks at some fraction of m max , and then decreases at higher slip velocities ( Fig. 1) (19)(20)(21)(22)(23). The decrease in  b /N, called rate-weakening drag, is caused by cavity growth with increasing slip velocity that shifts shrinking zones of ice-bed contact toward the tops of bumps. Owing to bump convexity, this shift reduces slopes of stoss surfaces in contact with ice. Laboratory experiments confirm this drag relationship for a sinusoidal bed (24).
The potential consequences of this slip law are severe: Weather or climate changes that increase slip velocity or decrease N would decrease resistance to slip, promoting faster flow and glacier instability, including glacier surge propagation (19,25). Despite convergence of process-based theories on this slip law, nearly all large-scale models of glacier motion parameterize drag to increase monotonically with slip velocity, with varying degrees of nonlinearity bounded by Coulomb slip in which  b /N is independent of slip velocity ( Fig. 1) (2,7,10,(13)(14)(15)(16). A major concern in applying the rate-weakening drag indicated by process-based theories is that most of them consider only 2D beds (19,20,22,23) with regular geometries. In the only 3D model, rateweakening drag was also found, but bed bumps were periodic and of uniform size and shape (21). Thus, the fundamental question of whether rate-weakening drag characterizes actual glacier beds is unanswered. Fig. 1. Slip laws for glaciers. Slip laws used in ice-sheet models (solid blue, dashed orange, and dotted purple) and a generalization of laws derived from process-based models of glacier slip that account for the effect of water-filled cavities at glacier beds (dashed light blue). orientation of preglacial discontinuities, such as fractures, bedding planes, and foliation (27,28). Thus, any effort to bracket the range of slip behavior for hard glacier beds needs to consider their geological and hence morphological variability.
On bedrock surfaces recently exposed by glacier recession, we measured diverse topographies of former glacier beds at high resolution with either a terrestrial laser scanner (TLS) or photogrammetry from an unmanned aerial vehicle (UAV) (see Materials and Methods). Proglacial areas selected for study, in Alberta, Canada, and in Switzerland (see figs. S1 and S2), consist of well-exposed sedimentary, metamorphic, and igneous bedrock with contrasting patterns of discontinuities. Associated bed morphologies range from meter-scale steps with planar stoss treads trending sub-perpendicular to the glacier-slip direction to bumps with convex stoss surfaces streamlined obliquely to the slip direction ( Fig. 2).
High-resolution surveys of proglacial bedrock span areas up to 0.7 km 2 , but because modeling glacier slip over 3D beds is computationally intensive, we develop sliding laws over smaller areas (140 to 270 m 2 ). These areas enclose subsections of the bed, which we call representative elementary areas (REAs), chosen to be morphologically representative of the full areas surveyed. REAs are identified using nine morphological characteristics that potentially affect sliding physics, and principal component analysis (PCA) is used to score the importance of each characteristic (see fig. S3). Scores for the full areas surveyed are compared to those for subsections of the bed to isolate REAs (Fig. 3). The REA for a particular surveyed area does not reflect its full morphological variability but does contain its most characteristic elements. Distinct morphologic differences, which largely reflect contrasting patterns of bedrock discontinuities, are evident among the REAs of the surveyed areas (Fig. 3).

Computed slip laws
To compute relationships between  b /N and slip velocity for the bed topographies of the REAs, we numerically model 3D flow of pure ice at its pressure-melting temperature in a boundary layer at the bed much thicker than the topographic relief (see Materials and Methods) (21). Ice at the top of the layer is driven at a uniform speed downglacier, as indicated by the average orientation of glacial striations on the bed, and cavities at the bed evolve to a steady size and shape (Fig. 4) under a steady, uniform effective pressure, N. Regelation is neglected. Owing to water that divides ice from rock, either in a thin film in zones of ice-bed contact or in cavities, local shear tractions on the sole of the ice layer are negligibly small. Therefore, the areally averaged drag at the bed,  b , is the integral of components of normal stresses on the ice sole parallel to the mean slip direction. Computing the stress distribution numerically is complicated by the irregular topography of the REAs, which rules out assigning at their edges periodic boundary conditions like those used in numerical simulations of glacier slip over simpler beds (20,21). To circumvent this problem, slices of bed topography measured directly adjacent to REAs are included in the modeling domain, with this additional topography tapered systematically near the domain edges to allow application of periodic boundary conditions (see Materials and Methods and fig. S4). We compute  b /N with a steady velocity applied to the  We determine the slip velocity, u b , in each case by integrating the velocity at the sole of the ice layer in the down-glacier direction.
Resultant slip laws are presented as scaled, dimensionless drag relations (19)(20)(21)23) of the form where n is the stress exponent, taken to be 3.0, in the power-law constitutive rule for ice deformation (29). The constant A s depends on the ice rheology and morphology of the bed and is defined on the basis of the long-standing drag relation for the case of no cavities (30) This scaling of slip velocity allows comparison among the REAs of differences in drag relations with and without cavities and results in a relation that is independent of whether the velocity or N is varied. Owing to the irregular morphologies of the REAs, which preclude scaling A s to a simple geometric attribute of the bed, the value of A s is computed through application of the model at velocities too small to induce cavity growth.
Conspicuously lacking from the derived drag relations (Fig. 5) is the significant rate-weakening drag ( Fig. 1) indicated by previous process-based models of hard-bedded glacier slip (19)(20)(21)(22)(23). For the REAs from the igneous and metamorphic beds surveyed,  b /N increases toward a bounding value with increasing scaled slip velocity (Fig. 5, A and B). For a sedimentary bed consisting of quasiperiodic steps with nearly planar stoss surfaces and steep lee surfaces, the bounding value of  b /N is reached at a very low slip velocity, so the drag relation closely approximates Coulomb behavior (Fig. 5C). The Coulomb behavior in this case is in contrast to mildly nonlinear relations between drag and slip velocity sometimes assumed for hard beds (7). A second sedimentary bed with more convex stoss surfaces indicates slight rate-weakening drag ( Fig. 5D) but also approximates Coulomb behavior. Repeating these calculations for these two sedimentary proglacial areas using different REAs yields drag relations with similar forms ( Forms of the derived drag relations are controlled by the slopes of stoss surfaces in contact with ice (Fig. 5). Although the maximum stoss slope, m max , controls drag for the case of idealized beds (18,23), for the irregular, 3D bed topography considered here, the maximum stoss slope in contact with ice is locally severe but spans an inconsequentially small area. We find that the average bed slope, parallel to the mean flow direction, in contact with ice scales predictably with  b /N as cavities get larger with increasing scaled slip velocity (Fig. 5). Rate-weakening drag, therefore, is not significant because stoss slopes of bumps in contact with ice tend to, with increasing slip velocity and cavity size, increase either slowly (Fig. 5, A and B) or rapidly (Fig. 5, C and D) toward limiting values rather than decrease. This effect can be due to stoss surfaces that are planar rather than convex (Figs. 2C and 3C). However, it arises more generally because actual 3D glacier beds are insufficiently periodic along a down-glacier flow line for cavities to commonly extend beyond the inflection points of convex stoss surfaces immediately downstream. This same lack of periodicity along a flow line also leaves some unusually steep stoss surfaces in contact with ice, which also inhibits rate-weakening drag (22).
Calculated magnitudes of  b /N are less than 0.2 (Fig. 5), which is less than steady-state friction coefficients for tills [0.32 to 0.84 (31,32)], but basal drag is generally greater for hard beds than for soft beds (31). This observation highlights the importance of subglacial hydrology. By impeding water drainage, soft beds promote higher subglacial water pressure than hard beds and thus have sufficiently low values of N to cause lower slip resistance than hard beds (17).

Hard-bed parameterization
These results indicate that a process-based slip model applied to observed glacier beds does not commonly produce rate-weakening drag. This finding qualitatively aligns the small-scale physics of hard-bedded slip with slip parameterizations of ice-sheet models that, on heuristic grounds, neglect the possibility of rate-weakening drag.
Our results also support the hypothesis that when water-filled cavities persist at the bed, the adverse slopes of parts of bumps in contact with ice control  b /N (Fig. 5) (18)(19)(20)(21)(22)(23)(24). This observation helps guide consideration of the effects on basal drag of bumps that are larger than the REAs and hence neglected in the modeling. Except near glacier margins where ice is thin, leeward cavities at scales greater than the REAs (i.e., tens of meters) will tend not to persist, owing to insufficient subglacial water discharge needed to sustain large, waterfilled, pressurized cavities. Nevertheless, on the up-glacier-facing slopes of larger bumps (i.e., "hills"), superimposed smaller bumps of the scales modeled herein will result in steeper up-glacier-facing slopes than in the absence of such hills. Maximum up-glacierfacing bed slopes set the value of Iken's bound (18) and C in Eq. 3 (20,21,23), so their values and values of  b /N would likely be somewhat larger than we have calculated (Fig. 5) if larger proglacial areas could have been modeled. However, because cavities at the scales of hills are expected to be absent or rare, the form of the calculated slip laws, which results from cavity growth with increasing slip velocity, is likely to be insensitive to hills on the bed.
The slip parameterization most commonly applied in ice-sheet models is Eq. 2, which is rooted in the physics of hard-bedded sliding without water-filled cavities at the bed (30). It provides a poor estimate of drag for the case of n = 3 (Fig. 5). Much larger values of n are sometimes used in ice-sheet models that use this parameterization (2,7,33), but it is then fully detached from process-based models of sliding physics because n effectively becomes a free parameter rather than having a value that, through the constitutive rule for ice deformation, depends on the properties of ice. Schoof (23), on the other hand, neglected the rate-weakening indicated by his 2D, process-based slip model to heuristically suggest a parameterization in which  b /N increases toward an upper bound.
To allow comparison with our calculated drag relations, we scale Schoof's parameterization [regularized Coulomb law (15)] as where C depends on the bed morphology (20). Although C cannot be readily estimated from the irregular topographies of the REAs, slip models require that C < m max (20,23). Using C as a fitting parameter with this constraint indicates that Eq. 3 can provide a good fit to our model results (Fig. 5).

Universal slip law?
Our results indicate that Eq. 3 is a reasonable parameterization for slip of glaciers over actual hard beds. Although the details of relevant processes are different, a slip law of the same form for soft beds (16,34) is supported by results of recent laboratory experiments on till (34). Moreover, application to Pine Island Glacier, Antarctica, of this slip law indicates that it accounts for decadal velocity changes substantially better than unbounded slip laws (i.e., Eq. 2) (15). Among applications of unbounded slip laws, those with sufficient nonlinearity to more closely approximate the regularized Coulomb behavior of Eq. 3 commonly best optimize fits of flow model results to remotesensing observations (33,35). Therefore, a slip law of the form of Eq. 3 may be universally applicable, allowing its parameters to be determined through application of observationally constrained glacier flow models without knowledge of the bed type (13). Given the sensitivity of modeled glacier discharges to various slip-law formulations, focusing on a single slip law and determination of its parameters would streamline efforts to estimate the contributions of glaciers to future sea-level rise.

Digital elevation model data
Digital elevation models (DEMs) for the four proglacial areas are constructed using data from a TLS and photogrammetry techniques with images collected from a UAV ( fig. S1). TLS data collected with a Riegl VA-400 scanner are used for the proglacial area of Castleguard Glacier. Scan positions are selected to assure scan overlap and mitigate shadow effects from bedrock features. The resulting point-cloud data are merged, filtered using reflectance (values ≤ − 25) and deviation (values ≤35) attributes, and manually cleared of spurious points with the RiSCAN Pro software package (36). For the three other proglacial areas (Schwarzburg, Rhône, and Tsanfleuron), photogrammetry techniques are applied to images collected from a UAV. The UAV is fitted with a 1-inch, 20-megapixel complementary metal oxide semiconductor camera that has an F2.8 wide-angle lens with a 9-mm focal length. The UAV adjusts its flight altitude to maintain a uniform distance from the ground while collecting images taken looking perpendicular to the ground surface. The images are used to generate point clouds with Agisoft photogrammetry software (37), which then are scaled and georeferenced with ground control points (GCPs) located across the surveyed areas. The 0.1, 0.05, and 0.1 m and total errors of 1.12, 0.03, 0.01, and 0.09 m for Schwarzburg, Rhône, Castleguard, and Tsanfleuron, respectively. Schwarzburg's total error is highest due to GPS failure partway through the photogrammetry survey.

REA selection
To determine an REA for each proglacial area, morphologic attributes of 10 m × 10 m subsections of a DEM are compared with those of the entire proglacial area. The DEMs are trimmed to a square and then subdivided into sections using a sliding window with a 2-m step spacing. The subsections are detrended, and their morphologies are characterized after Leach (38) and Stout et al. (39) with the following nine parameters: root mean square (RMS) roughness, RMS slope, and RMS curvature to characterize the degree of surface roughness; texture aspect ratio to characterize the isotropy of the topography; length asymmetry and slope asymmetry to characterize the asymmetry of bumps in the direction of ice flow; skewness to characterize the symmetry of the surface elevation distribution; kurtosis to characterize the sharpness of bumps; and a 2D correlation coefficient, r, between the Fourier spectrum of the proglacial area and that of each subsection. A description of how these parameters are calculated is provided in the Supplementary Materials (text S1).
To find the REA of a proglacial area, PCA (40) is performed on the nine parameters to help identify a subsection that is representative of the DEM (see the Supplementary Materials for more details about PCA). In contrast to multidimensional distance methods (e.g., Euclidean, cosine, or Mahalanobis) that give equal weight to all parameters regardless of the variance they capture, the PCA method determines the parameter weights to capture the most variance of the proglacial morphology dataset. As the parameters that best describe the DEM morphologies are unknown beforehand, this approach allows the determination of the most effective variables (PCs) for describing the morphology and the percentage of the total variance captured by these variables. To evaluate the level of similarity between a subsection and the DEM, a representative value for the entire proglacial area, X  , is designated for each parameter, where  is the parameter index. For r, X  is set to one, which represents a perfect power spectrum match between a subsection and the entire DEM (see text S1). For all other parameters, X  is set to the mean of the parameter set. The X  values are normalized with the other parameter values to a mean of zero and variance of one before performing the PCA, owing to the different units of parameters. The PC scores (magnitude of each PC in the PC coordinate system) are estimated for the normalized parameter values of individual subsections,  ,i , and the representative parameter values of the entire proglacial area,   , where i denotes the subsection index. The representative PC scores,   , are subtracted from each subsection's score,  ,i , resulting in a similarity index,  i , for each subsection that is used to find the most representative subsection of the DEM Here, values closer to zero reflect a greater degree of similarity between a given subsection and the entire proglacial area. Owing to the use of PCs, this similarity index intrinsically incorporates the appropriate weights of each parameter to best describe the surface attributes. Following Cattell and Vogelmann (41), the number of PCs used to locate an REA is determined by finding the point at which the percentage of the total variance captured by successive PCs converges to an approximately constant value and then including one PC beyond that point. For all of the proglacial areas, the subsections with the lowest combined  value from the first two PCs are taken to best represent the morphology of the proglacial areas surveyed, as the percentages of the total variance of components beyond the first PC are relatively uniform. Loadings (level of correlation of each parameter with a PC) of the first two PCs for the four different DEMs are shown in fig. S3. The first two PCs account for 41, 42, 46, and 43% of the total variance for the proglacial areas of Schwarzburg, Rhône, Castleguard, and Tsanfleuron, respectively.
Application of periodic boundary conditions in the numerical model requires choosing a subsection with minimal elevation differences between its opposing sides. For each proglacial area, we consider the 50 subsections (∼0.5% of the total number of subsections of a proglacial area) with combined similarity indices closest to zero (greatest similarity) as a set of potential REAs. Each subsection of this set is expanded by 5 m in each direction using the surrounding bedrock topography, resulting in a 20 m × 20 m section centered on the statistically representative 10 m × 10 m subsection. As a measure of periodicity of any given rectangular section, we consider the L 2 -norm of the differences in elevation between opposing sides averaged over the length of each boundary. In an expanded section, this periodicity measure is calculated for every possible rectangular section centered on and containing the subsection. This procedure is performed for each expanded section in the set of potential REAs, with the rectangular section with the highest periodicity chosen as the section best suited for the numerical model. In cases where surface debris obscures bedrock, the rectangular section with the next highest periodicity is used. Last, a low-pass Butterworth filter of order six is applied to the selected rectangular section to remove wavelengths smaller than 0.25 m, making the shortest wavelengths present in the topographical data similar to the horizontal resolution of the numerical model.
To make the boundaries of a selected rectangular section fully periodic, we use a periodicizing tapering algorithm with a taper width of 1 m at the boundaries (see text S2). A 2D rectangular sliding window is used to average the elevation values within the tapered zone. The maximum size of the window is 2 m in each direction, with the window centered on the point (elevation value) that is averaged. The size of the window progressively decreases to zero as averaged coordinate points approach the inner edge of the tapered zone (see fig. S8). The varying size of the averaging window results in a smooth transition from the tapered part of the section to its untapered interior. In choosing a modeling domain, for the periodicity measure, we considered rectangular sections only large enough so that the centered REA is not affected by the tapering algorithm. For the four proglacial areas, comparing the PC scores of the tapered and low-pass-filtered rectangular section containing the selected REA to the corresponding unmodified section indicates that only a small statistical change in bed morphology is introduced by the tapering procedure. The relative change in the PC scores for the selected rectangular sections of the proglacial areas of Schwarzburg, Rhône, Castleguard, and Tsanfleuron are 0.07, 16.44, 1.31, and 0.98%, respectively. The high relative change of the Rhône proglacial area is due to the long wavelengths of its bed morphology and the small area of the selected rectangular section.

Modeling strategy
The equations governing the viscous deformation of ice and steadystate geometry of the glacier sole, z s , are given by the Stokes equations (5) and kinematic free-surface equation (7). For the model, ice occupies a Cartesian domain in a boundary layer at the glacier bed, which consists of the rectangular domain that includes the REA and surrounding tapered regions. To the top of the domain, we apply a flow velocity, u e , and ice overburden pressure, p i . Gravity in the thin boundary layer can be neglected. At the bottom of the domain, the glacier sole either is in contact with the underlying bedrock, z b , or is the roof of a water-filled cavity, z c . The bedrock is taken to be impenetrable and the water pressure in cavities, p w , to be spatially uniform over the bed. Water in a thin film divides ice at its pressure-melting temperature from rock where water-filled cavities are absent, so the low viscosity of water results in negligible tangential stresses everywhere at the glacier sole (free slip). Basal melting and regelation are neglected. At the vertical boundaries of the tapered domain, periodic boundary conditions apply, so ice deforms during slip over morphologically repetitive bedrock of infinite extent.
For a particular bed and cavity geometry, the velocity, u = (u x , u y , u z ), and pressure, p, are then given as the solution to ∇ ·u = 0 (conservation of mass) (5b) n ·  · n = − p i , u x = u e , u y = 0 (top boundary condition) (5c) n ·  · t = 0 on z s (free slip condition) (5d) u · n = 0 on z b (impenetrability condition) (5e) n ·  · n = − p w on z c (cavity roof condition) (5f) where  and  are the Cauchy and deviatoric stress tensors, and n and t are the outward-pointing normal and tangent to the boundary, respectively. The deviatoric stresses are related to strain rates through the constitutive equation for ice deformation (29) where  is the effective ice viscosity dependent on strain rate, with n = 3. In general, the rate factor B depends on temperature but is constant in the isothermal case that is appropriate here (42,31). In Eq. 6, the strain-rate tensor and effective strain rate are given by  ̇ ij = 1 _ 2 (∂ u i / ∂ x j + ∂ u j / ∂ x i ) and  ̇ E = 1 _ √ _ 2 (  ̇ ij  ̇ ij ) 1/2 , respectively. The Stokes equations (5) describe steady flow, but steady-state geometries of cavities at the bed are not known a priori. Thus, a separate equation is needed to simulate the transient evolution of z s to the desired steady state in which water-filled cavities at the glacier sole are not shrinking or expanding. Considering the glacier sole to be a free surface, with the additional condition that it be at or above the underlying bed at every point, the kinematic free-surface equation is ∂ t z s + u x ∂ x z s + u y ∂ y z s = u z , such that z s ≥ z b For a steady-state ice geometry,  b is the bed-averaged viscous drag in the general flow direction (19,22,23), calculated as where  is the area of the computational domain in the horizontal plane and n x is the x component of n. In a corresponding manner, the average slip velocity in the direction of flow is calculated as The set of Eqs. 5 and 7 is solved numerically using the finite element software Elmer/Ice (43). All simulations are performed using nodal finite elements on meshes generated by extruding an unstructured triangulated (x/y plane) base mesh in 30 layers in the vertical (z) direction, resulting in prismatic/wedge elements. The same vertical mesh extrusion is used for all glacier topographies, with the density of element layers highest near the base of the computational domain. The top of the domain is fixed at z = 15 m, whereas the bottom of the domain adjusts to the geometry of z s in each time step. All base meshes are generated with a characteristic cell size (triangle edge length) of 0.15 m using the frontal Delaunay algorithm through the software Gmsh (44). Owing to the different horizontal extents of the computational domains, this results in meshes consisting of 938,640 (Castleguard), 518,520 (Rhône), 962,460 (Schwarzburg), and 777,840 (Tsanfleuron) cells. Additional description of the numerical model is given by Helanow et al. (21) and its supporting information, where results of sensitivity and parameter tests can also be found.