Research ArticleGEOLOGY

Sliding dominates slow-flowing margin regions, Greenland Ice Sheet

See allHide authors and affiliations

Science Advances  10 Jul 2019:
Vol. 5, no. 7, eaaw5406
DOI: 10.1126/sciadv.aaw5406


On the Greenland Ice Sheet (GrIS), ice flow due to deformation and sliding across the bed delivers ice to lower-elevation marginal regions where it can melt. We measured the two mechanisms of motion using a three-dimensional array of 212 tilt sensors installed within a network of boreholes drilled to the bed in the ablation zone of GrIS. Unexpectedly, sliding completely dominates ice motion all winter, despite a hard bedrock substrate and no concurrent surface meltwater forcing. Modeling constrained by detailed tilt observations made along the basal interface suggests that the high sliding is due to a slippery bed, where sparsely spaced bedrock bumps provide the limited resistance to sliding. The conditions at the site are characterized as typical of ice sheet margins; thus, most ice flow near the margins of GrIS is mainly from sliding, and marginal ice fluxes are near their theoretical maximum for observed surface speeds.


Sliding and deformation of the Greenland Ice Sheet (GrIS) exert a first-order control on ice sheet mass balance by dictating ice transport to marine outlets and low elevations prone to high melt rates. However, while the surface velocity of the GrIS is well documented at dense spatial (1, 2) and temporal scales (38), the relative contribution of the two components of motion remains uncertain and is known in a few locations only (911). Estimates of sliding speed are typically based on the residual between observed surface velocity and modeled ice deformational velocity (1214), but uncertainties related to the stress field and the constitutive law for ice inhibit high confidence in the basal velocity derived from data-constrained model inversions (15). Precise knowledge of ice temperature is rare, and evidence exists for chemical impurities (16, 17), variable grain size distributions (1618), and preferred ice fabric (17, 18), all of which have potential to strongly affect the constitutive law for ice. This creates large uncertainty in partitioning the contributions of sliding/deformation to surface velocity: past estimates of ice flow across the GrIS ablation zone contrast, having suggested that motion is controlled either by sliding (12, 13) or by deformation (14).

Although surface speeds and presumably basal motion are far lower during winter than the melt season, winter flow conditions persist through the majority of the year and contribute the most to the net annual displacement on the GrIS margin (6, 19). Consequently, winter ice flow, despite the lack of surface meltwater forcing, is the central driver of the annual mass flux along the margins of the GrIS (6, 19). Much of the uncertainty in the interpretation of annual surface speeds and mass flux thus stems from the fact that confident partitioning of winter sliding and deformation is limited to very few measurements (9, 10).

We address this fundamental problem with a unique and comprehensive set of in situ measurements and analysis of ice motion. Our data are recorded by a network of GPS stations and a dense three-dimensional array of englacial inclinometers installed in a cluster of boreholes drilled in the ablation zone of the GrIS. Our measurements are used to explicitly partition the surface velocity into its deformation and basal sliding components during the 2015/2016 winter (~240-day period between melt seasons) and to elucidate the processes controlling winter motion along the basal boundary. Unaccounted for in previous inclinometer experiments on the GrIS (911), the density and scale of our multidimensional array overcome possible biases that arise from localized deformation variations (20). This yields a high-fidelity representation of deformation and sliding velocities over a length scale of an ice thickness.

The field site is located 33 km from the terminus of Issunguata Sermia within the ablation zone of the western margin of the GrIS (Fig. 1). The region is land terminating and is characterized by surface speeds on the order of 100 m a−1 (Fig. 1). Our instrumented block of ice is ~3.5 × 107 m3 and occupies a topographic plateau (Fig. 1 and fig. S1) (21). The ice is 640 to 688 m thick, is temperate at its base, and reaches a minimum of about −12°C mid-way through the ice depth (22). The borehole grid horizontal span is nearly an ice thickness, with dimensions of 465 m along the flow direction and 280 m across flow, and is bounded by a five-station GPS network. A variety of experiments conducted in the boreholes indicate that the bed is hard bedrock with no thick and deforming sediment layer (23).

Fig. 1 Study area and field site setting.

(A) Three-dimensional representation of surface topography, bed topography (21), borehole locations (vertical colored lines), and GPS locations (black surface markers). (B) Regional surface velocity map showing winter flow conditions in 2010 (2, 43). Gray contours indicate bed topography in meters above sea level (m a.s.l.) (21). Left inset shows location of study region. (C) Surface map of study location showing GPS stations (black triangles, text—mean winter velocities) with flow vectors and borehole locations (colored circles: 14W, purple; 14N, blue; 14Sa, green; 14Sb, red; 15Ca, gray; 15Cb, cyan; 15S, pink; 15E, orange; 15N, dark gray). Surface strain rate tensor calculated across the GPS network: ε˙xx = 0.0006 a−1, ε˙yy = 0.0037 a−1, ε˙xy = 0.0014 a−1 .

Measured deformation and sliding

Deformation. Ice deformation within the study block is dominated by along-glacier shear and is limited to ice ~150 m above the bed (mab); almost no deformation occurs in the upper 75% of the ice column (Fig. 2). The shearing region corresponds to warmer ice (i.e., ~0° to −7°C) near the bed, although the temperate basal layer occupies only the lowermost few meters. The vertical profile of deformation is strongly variable from place to place around the block (Fig. 2). Some profiles show almost no deformation from surface to bed, and other areas reach maximum shear deformation rates of up ~0.15 per year. The largest deformation rates are not observed closest to the bed, but rather peak values are several tens of meters above the bed, thus forming a characteristic “kink” in the vertical deformation profile (described in more detail below).

Fig. 2 Englacial measurements of temperature and ice flow.

(A) Englacial temperature profiles measured at each borehole location. (B) Site mean deformation profile (black markers), plotted with strain rate measurements in background. Pink envelope encompasses all the shear strain rate observations. Calculated deformation boundaries (gray envelope) given in situ ice temperatures and an estimated range of driving stress (76 to 107 kPa). (C) Velocity profile (blackline) derived by integrating site mean deformation profile. The surface velocity (us) is partitioned into basal sliding (ub) and ice deformation (ud) components.

A calculated Glen-type deformation field for the block, representing a range of driving stresses, has limited overlap with our observations (Fig. 2). The driving stresses are computed at greater than ice thickness length scale (24) as 76 to 107 kPa (Materials and Methods) in the direct vicinity of the site—values that are quite typical of glaciers and ice sheets (24) and are representative of large areas of western Greenland (25). Further, our temperature field across the block (Fig. 2) (22) is used to constrain the viscosity parameter’s strong temperature dependency. The individual measurements and the mean deformation profile of the entire ~3.5 × 107 m3 block of ice are not described by the field of Glen-type deformation, with standard flow law parameters. The calculated along-glacier shear extends about twice as high upward into the ice as observed, the calculated magnitudes are up to twofold too large, and our data show irregular deformation profiles. The discrepancies could be due to deformation mechanisms poorly represented by Glen’s law, properties of the ice such as fabric making it stiffer than expected, and/or unaccounted-for higher-order stresses. The presence of extensional normal strain rates at the surface (Materials and Methods) suggests some contribution from higher-order stresses to the stress state of our instrumented block. However, the fact that the driving stress we calculate is representative of western Greenland (Discussion) indicates that the stress conditions are not anomalous and the influence of higher-order stresses is small.

Sliding. During the 2015/2016 winter, the mean surface velocity was 114.8 m a−1 (±1.3 × 10−6 m a−1 SD). This rate is typical of land-terminating regions of western GrIS and is an order of magnitude less than some fast-moving outlet glaciers (2). Our inclinometry data permit direct partitioning of the surface velocity into mean ice deformation and basal sliding components (Materials and Methods). We find that ice deformation contributes 4.6 m a−1 (0.4 to 10.8 m a−1 uncertainty bounds) and basal sliding is responsible for 110.2 m a−1 (104.0 to 114.4 m a−1 uncertainty bounds) (Fig. 2). Thus, only 4% of the surface motion can be attributed to deformation of the ice column, with basal sliding producing the other 96%. This is best described as plug flow, where the surface motion is overwhelmingly derived from basal motion.

Basal boundary. The temporal resolution and spatial density of ice deformation measured by inclinometer tilting near the bed permit detailed investigation of ice motion along the basal boundary. As the sensors moved approximately 70 m across the bed during the winter period, measurements of inclinometer tilting demonstrate variations: variations within individual records, and variations between different inclinometers located at similar ice depths. These could potentially arise from processes varying in time, space, or both (Fig. 3 and fig. S3 and S4). However, the variations between sensors in the ice tens of meters above the bed are uncorrelated between borehole locations or to the surface velocity measured during the winter (fig. S4). While transients cannot be entirely discounted (Materials and Methods), this suggests that the variations in the tilt rate are not driven by widespread processes varying in time, where a response in ice speed would be expected, but instead are best attributable to localized spatial variations in the basal stress field recorded as the inclinometers slide across the bed.

Fig. 3 Observed and modeled basal tilt rates.

The tilt rates (5 day) for the basal inclinometers in 15Cb (A), 15S (B), and 14Sa (C). Top panel for each hole is the tilt rates for the bottom five inclinometers plotted versus basal distance traveled during the 2015–2016 winter. The line color corresponds to the inclinometer installation height in meters above the bed (mab), as plotted in the bottom panel. The bottom panel shows the same data plotted as interpolated surface. The right column shows the best-match modeled tilt rates (D to F) for inclinometers sliding across beds with different roughness characteristics. The model type and roughness are indicated in the plot title. The top panel shows modeled inclinometer tilt rates plotted against basal distance, and the bottom panel shows the magnitude of the modeled tilt rates adjacent to the bed and the flow paths (colored lines) for the modeled basal inclinometers. The line colors in the top panel correspond to the inclinometer flow paths in the bottom panel. Each model scenario plotted is the best match to the adjacent observed tilt rate patterns. Model outputs for all model scenarios are plotted in fig. S5.

Transforming the tilt records to their spatial dimension reveals the pattern of ice deformation as the study block slides over the bed (Fig. 3 and fig. S4). Directly adjacent to the bed, where the highest shear deformation is expected, unusually low and reverse tilting is measured by the bottom inclinometer in some sensor strings (14Sa, 15Cb, and 15S). Records from these strings demonstrate sudden transition from negative tilting (up-glacier shearing) to positive tilting (down-glacier shearing) as the sensors slid over regions of the bed. These variations occur over 5- to 40-m horizontal length scales (Fig. 3). Tilting variations are also observed in sensors extending 15 to 45 mab, with diminishing amplitude moving upward. The spatial patterns of deformation higher in the ice column are nearly out of phase, offset, or uncorrelated to the tilt rate patterns measured closest to the bed (Fig. 3). Overall, these variations reveal large spatial changes in the basal shear stress at the tens of meters scale. Averaged over the study block, these spatial variations result in the pronounced basal kink in the mean vertical profile of along-glacier shearing (Fig. 2).

Analysis of sliding

Sliding over a rough bed. On a hard bed, the spatial heterogeneity in tilt rates that we measure can occur because of bed roughness and/or spatial changes (i.e., “patchiness”) in basal slipperiness. A suite of simple flowline modeling experiments (Materials and Methods and table S2) was performed to investigate whether sliding over a range of simplified basal boundaries could reproduce the spatial variations observed by the near-bed inclinometers. The basal boundaries are a synthetic representation of different types of bedrock roughness, patchiness (alternating sticky and slippery patches), and combination roughness/patchy beds (fig. S5 and table S2).

Model results show that the inclusion of bed roughness is required to generate the negative tilt rates and spatial tilting patterns, consistent with our measurements near the bed (Fig. 3). At borehole 15Cb, sinusoidal roughness provides the best match to the observed variations. Sliding over isolated bedrock asperities most closely reproduces measurements at holes 14Sa and 15S. In these scenarios, modeled low–tilt rate regions occur locally within topographic lows preceding or following the roughness feature, and high tilt rates occur above the roughness peaks (Fig. 3). Inclusion of patchy slipperiness with bedrock roughness does not appreciably change the modeled pattern of tilt rates from the model scenario prescribing bed roughness alone (fig. S5), indicating that variations in slipperiness have a second-order influence on deformation variability near the bed. The magnitude of the tilting variations is higher in the modeled scenarios than observed. This could reflect that either driving stress (100 kPa) or rheology used in the modeling scenarios is either higher or softer than that of the field location. Overall, the results show that the measured variations in the near-basal tilt rates at the field location are most consistent with rapid sliding over bedrock roughness with characteristic length scale on the order of tens of meters.

Following this interpretation, high tilt rates over bed asperities indicate that a substantial portion of the basal resistance is accommodated by sparsely spaced topographic features tens of meters in scale. Away from these features, large areas with negative tilt rates suggest that much of the bed supports little resistance to sliding. In modeling scenarios, this manifests in parameterized slipperiness, which, to achieve similar flow speeds as the study block, must be tuned so that the ice-bed interface supports little shear stress, and stresses build up on the stoss side of roughness features (Fig. 3, fig. S7, and table S2). This is unexpected as smaller wavelength features (~1 m), which are expected to take up much of the resistance at the bed (26), are noticeably absent. It is possible that deformation effects from small-scale roughness are aliased by our measurements because basal tilt sensors are not immediately above the bed. However, in hole 14Sa, the basal sensor is installed just 2 mab, and measured tilt variations only occur from what we interpret to be a bed asperity that is 10 to 20 m wide. However, given that meter-scale roughness is expected (23), we posit that limited ice-bed contact from substantial cavitation could significantly lessen the basal resistance of small-scale features. This provides a reasonable mechanism for both the high sliding fraction and the concentration of stress on large bed asperities as observed.

Integrated effect on deformation. Averaged over the wavelength of the bed roughness, modeling results show that rapid sliding over a rough bed generates a shear strain rate profile through the ice column that does not smoothly increase with depth, but instead has a sharp kink in basal ice—as observed in our instrumented block. This finding persists in all bedrock roughness scenarios (fig. S6). The modeled deformation kink arises due to the vertical redistribution of shear stress by higher-order stresses, which act to enhance shear stress, and therefore shear strain rates above bed asperities, and resist flow in topographic lows between bed roughness features (fig. S7). Sliding over a slippery and rough bed is therefore self-consistent with the following various aspects of the data: the rapid sliding, detailed tilt sensor measurements, and deformation profiles we measure at the site.

Our analysis thus implies that basal kinks are characteristic of sliding over a topographically complex basal boundary at a length scale of tens of meters. Similar deformation kinks directly above the basal interface have also been observed in the few existing in situ measurements of Greenland ice deformation (911), including over till beds (9, 11). This suggests that at least some fraction of the basal motion over till beds is resisted by bed roughness, demonstrating the ubiquitous importance of hard-bedded sliding physics to Greenland ice flow.


Rheologic factors could vary widely (14, 16, 17, 9, 10, 27, 28) in western GrIS, and in the absence of in situ constraints, calculations of deformation velocities can span orders of magnitude (Supplementary Materials and Methods and fig. S8). While ice in western GrIS has been suggested to be predominantly Holocene in age (29), evidence from locations toward the margin from our site indicates that pre-Holocene ice is present (30) and this ice should be prone to enhanced deformation (16). This, and the absence of weak basal sediments, may imply substantial winter displacement from ice deformation: the opposite of the plug flow we observe. Our measurements thus show that pre-Holocene ice is either too thin or not soft enough to have substantial impact on ice flow given the stress conditions at the site. Perhaps closer to the margin, where temperate ice occupies a larger fraction of the ice column (27), more of the driving stress is accommodated by deformation of deep warm ice and sliding is less prevalent. However, most of the ablation zone lies inland of our site, with even colder ice (27) and a comparable driving stress (Fig. 4) (25).

Fig. 4 Driving stress across the western margin of the GrIS.

Driving stress calculated for the ablation zone of the western margin of the GrIS. Dots indicate field locations with direct measurements of deformation and sliding. Red dot is our field location, and cyan dots show other field locations with direct partitioning of surface motion.

The few previous measurements of sliding/deformation partitioning were made over beds with thick sediment cover (31, 32) and/or in regions strongly influenced by fast-flowing marine outlet glaciers (9, 11), where a strong proportion of basal sliding is expected. A sliding fraction of 44 to 73% was measured during the winter (10) over a bed identified as a thick layer of sediment (31). On marine-terminating Store Glacier, which is underlain by weak sediment (32), sliding accounted for 63 to 71% of the summer surface speed (11). In addition, along the flanks of Jakobshavn Isbrae, a basal sliding fraction of 63% was measured (9). Yet, we find ice flow that is almost entirely by basal sliding in arguably the least likely setting to show high winter sliding. Further, our measured sliding conditions are not anomalous to a single profile but are characteristic of a three-dimensional ice volume. Direct observations of basal sliding now span the spectrum of anticipated physical conditions at the bed. Only sliding-dominated flow has been observed, including during the winter over a hard bed.

Marginal environments are also topographically complex, including both ridges and troughs (21). Deformation and sliding measurements now span these topographic settings, indicating that sliding-dominated motion occurs in troughs with relatively thick ice (10, 11) and also over topographic highs where thinner ice is present (9).

Sites where deformation has been measured in western GrIS not only span the physical conditions of the bed and topography but also cover a broad range of driving stress and winter surface velocity conditions (Figs. 4 and 5). In this context, we find that our measurement of sliding-dominated ice flow in winter is representative of typical marginal ice flow conditions. Given the aggregated evidence, we hypothesize that annual surface displacement mainly represents sliding motion, and deformation plays a limited role in flow on the land-terminating margins of the GrIS.

Fig. 5 Driving stress versus surface velocity, western GrIS.

Driving stress calculated across the western margin of the GrIS plotted against collocated, InSAR-derived surface velocities collected between 1 December 2008 and 28 February 2009 (2). Gray dots show all the measurements, while the red circle marker is measured from the field location and the blue dots represent those measured from the other locations with direct partitioning of surface motion (17, 18, 26).

Hard-bedded regions serve as the conceptual basis for our understanding of sliding and drainage physics on the GrIS (33, 34) and are representative of select regions of the GrIS (23, 35). Our results show that a plug flow–like behavior exists over a hard bed during the winter months, suggesting the paradigm of hard-bedded regions as typically low sliding should be reconsidered. Moreover, our analysis demonstrates the importance of the bed roughness at the tens of meters scale for influencing ice deformation and affecting overall basal coupling. That we find winter ice flow to be almost entirely by sliding over a hard bed and outside of fast-flowing regions of GrIS is unexpected, but we show that our site is representative of the margins in terms of driving stress and surface velocity. With all direct evidence now pointing to sliding-dominated flow across GrIS’s ablation zone, ice flux along the flanks of the ice sheet should be considered near its theoretical maximum for observed surface speeds. This suggests that the GrIS can efficiently redistribute mass and move ice to the low-elevation and high–melt margin zones and respond more rapidly to changes in climatic boundary conditions than previously known.


Inclinometer construction and installation

Individual inclinometer probes were constructed using two orthogonal high-precision analog inclinometers (Murata SCA103T) that measure tilt from vertical (θ). Each inclinometer has a resolution of 5 × 10−4° and a range of ±15° from vertical. The analog output from the inclinometers is measured by a separate 24-bit A/D converter. The orientation of the inclinometer is measured by a three-axis magnetometer (Honeywell HMC5883) that has a resolution of 2°, and is laboratory calibrated with an eight-point calibration to ±10°. Given the increased uncertainty of the magnetometer near true polar north, the field resolution of the magnetometers is estimated at ±20°. The inclinometer probes also have an onboard temperature sensor to measure ice temperature. The temperature sensor (Texas Instruments TM102) has a resolution of 0.0625°C, is factory calibrated to ±1°C, and is laboratory- and field-calibrated to ±0.25°C. Each inclinometer package is embedded within a cylinder using epoxy and constructed into inclinometry strings using CAT-5 cable.

A total of nine boreholes were drilled from the surface to the bed at the field location (Fig. 1) using hot water methods (36) during the 2014 and 2015 melt season. The boreholes were drilled in a grid, with an E-W and N-S extent of 465 and 280 m, respectively (Fig. 1). The borehole locations are denoted by the year they were drilled and the relative location within the borehole grid (east, west, north, south, and center). Strings of inclinometer probes were installed in each borehole immediately after drilling. The strings consisted of inclinometer probes spaced 10 m apart from 5 to 145 mab and 20 m apart from 145 to 625 mab. At two borehole locations, 14Sb and 15Cb, only the bottom half of the ice column (<165 and <325 mab, respectively) was instrumented with inclinometers. After installation, refreezing of the borehole coupled the inclinometer probes to the ice. No thick temperate layer was measured, with ice temperatures recorded below the pressure melting point by the basal most inclinometer at each string. Once installed, data were recorded continuously from each inclinometer string using a custom surface logger at 2- or 4-hour intervals. One borehole, 15N, yielded no usable inclinometry data after installation.

Winter deformation rates

Continuous tilt measurements (θ) recorded by each inclinometer probe were used to estimate the vertical gradients in horizontal velocity (du/dz) across the 2015–2016 winter. The following relation estimates the shear strain rate from inclinometer tilt in the x-z planedudz=ΔxΔzΔtΔtanθaΔt(1)where Δt spans the length of the winter period, and θa is tilt in the direction of ice flow.

The tilt of an inclinometer installed with its cylindrical axis oriented vertically is solely due to vertical shearing of the ice column (37). Once rotated from this vertical orientation, analytical modeling of inclinometers embedded in deforming ice has shown that the rotation rate is influenced by normal strains (37). This can affect the estimate of du/dz using Eq. 1 if field conditions deviate significantly from shear dominated flow. The magnitude of the normal components of the strain rate tensor calculated from the GPS surface array (ε˙xx = 0.0006 a−1, ε˙yy = 0.0037 a−1, ε˙zz = −0.0043 a−1, and ε˙xy = 0.0014 a−1) shows that ice deformation is not only due to vertical shearing. This suggests that it is necessary to quantify the impact of normal strains on our estimates of du/dz using Eq. 1.

To do so, we used a sensitivity analysis (Supplementary Analysis) that shows that the tilt rate of near-vertical inclinometers is dominated by the rate of vertical shearing. Specifically, the analysis illustrates that when θ < 3°, Eq. 1 is a robust method to estimate du/dz. When θ < 10°, Eq. 1 is still reliable for estimating du/dz under most conditions expected at the field location. As a metric of data quality, we therefore partition deformation measurements into two sets based on whether the mean inclinometer tilt over the winter period is θ < 3° or θ > 3° when presented in fig. S2.

The instrumental error in the calculated vertical shear strain rates ranges from ±0.5 × 10−5 to 1.7 × 10−5 a−1 (SD). Inclinometers with nonfunctional tilt sensors or those with sensor errors identified by the logger during sampling were not included in the analysis. Uncertainty associated with the effect of normal strains is estimated to be an additional ±10% of the strain rate estimated using Eq. 1 (Supplementary Analysis).

The winter period spans 240 days during the 2015–2016 winter [day of year (DOY) 2015, 255 to 495]. During this period, the velocity field did not undergo any large and rapid variations in time (fig. S4). Air temperatures generally remained below freezing, with the exception of a brief excursion to positive temperatures near the end of the winter (DOY 2015, 466 to 471), which resulted in a minor velocity decrease at the westernmost GPS station. The lack of coherence between GPS stations and small magnitude of the event suggests that this is likely due to coupling to nonlocal flow downgradient of the site.

For most inclinometers, du/dz was estimated using Δt equal to the entire winter period. Discontinuous or intermittent records at specific borehole locations and inclinometers yield some variation in the length of the record used to estimate du/dz. This is documented in table S1.

Deformation and basal sliding velocity

The deformation velocity, the component of motion from shearing of the ice column, was calculated by integrating the du/dz estimated from the inclinometers from the surface to the bed. du/dz values at each inclinometer depth were averaged across all boreholes to account for any local variability present in the deformation rates. The result is a profile of mean deformation rates that is representative of the field location. There is a small offset in the inclinometer installation depth (2 to 4 m) in holes 14Sa, 14Sb, 14W, and 14N compared to the inclinometer strings installed in 2015. In the averaging process, du/dz values in boreholes drilled in 2014 were averaged with du/dz values from the closest depth interval from boreholes drilled in 2015.

Integration of the mean deformation profile yields an estimate of the deformation velocity for the field site. du/dz at the bed was set to the value calculated from the bottom inclinometer depth, and du/dz at the surface was set to zero. Measurement gaps between inclinometer depths were filled using linear interpolation. Uncertainty bounds in the deformation velocity were estimated by calculating the deformation velocity using the minimum and the maximum du/dz measured at each depth interval at any borehole location. This yields a conservative range for the error in the deformation velocity beyond that of the averaging and integrating the individual sensor error and uncertainty.

The basal sliding velocity was calculated as the residual of surface velocity and deformation velocity. The surface velocity was calculated from position data collected continuously at 15-s intervals from five Trimble Net-R9 GPS receivers installed in a diamond array (Fig. 1). The data from the GPS receiver were processed against a base station located ~22 km away off the ice sheet using TRACK version 1.29 differential kinematic processing software. The winter surface velocity at each GPS location was calculated usingus=ΔxΔt(2)where Δt is the duration of the winter period and Δx is the change in position across the winter period. The error in the velocity is estimated to be ±1.3 × 10−6 m a−1 (SD) for the velocity at each station. The site mean velocity was calculated by averaging the mean winter velocity from all five GPS locations. For continuous velocities presented in fig. S4, the position data were smoothed using a Gaussian smoothing kernel, with a window length of 5 days, and then differentiated continuously through the winter period with Δt = 5 days.

Ice thickness and driving stress

The driving stress was estimated using field measurements of ice thickness and surface slope calculated from a surface digital elevation model (38). Measurements of ice thickness were made at each borehole location during drilling using marked drill hose and instrument cables. Measured ice thickness ranged from 640 to 688 m across the field site, with a mean ice thickness of 668 m. The elevation data were smoothed using a Gaussian smoothing kernel over an area of 1 km2. The surface slope was calculated across a distance of 1 km, yielding values ranging from 0.74 to 1.04° in the direct vicinity of the field site. The driving stress was then calculated using the followingτd=ρigHtan(α)(3)where ρi is the ice density (900 kg/m3), g is the gravitational constant (9.81 m/s2), H is the mean ice thickness, and α is the range of the surface slopes for the field site.

Calculated deformation rates

Theoretical deformation profiles were generated using in situ measurements of ice temperature and the driving stress estimated for the field location (Fig. 2). Deformation rates through the ice column were calculated using a Glen’s style ice rheology given bydudz=2A(T)(τddH)n(4)where A(T) is a temperature-dependent creep parameter, d is the depth below the surface, and n is the flow exponent. The above formulation assumes a simplified flow field where deformation rates are solely a function of the driving stress. Using recommended values (24), the temperature-dependent creep parameter is calculated through the full thickness of the ice column using the temperature profile collected at 15Ca. The flow exponent is assumed to be n = 3, the most commonly accepted value used in glaciology (24). No basal enhancement factor was applied in the calculation of theoretical deformation rates.

Continuous tilt rates versus distance

Continuous tilt rates for inclinometers less than 55 mab were calculated for 14Sa, 14N, 15Ca, 15Cb, 15S, and 15E. Inclinometry strings installed in holes 14Sb and 14W had discontinuous winter data, and continuous tilt rates were not calculated. The tilt time series from each inclinometer axis was preprocessed with a Hampel median filter, with a window length of 2 days to remove spurious signal noise and outliers. The tilt measurements were then smoothed using a Gaussian smoothing kernel with a window length of 5 days. Tilt rates of the inclinometers were found by calculating the change in tilt of the sensor with respect to time using the following relation, where Δt = 5 daysθ˙=θ1θ0Δt(5)

The tilt rates are plotted versus distance traveled relative to the bed in Fig. 3 and fig. S3. This distance was estimated using the winter basal sliding velocity to convert each inclinometer time series to a distance domain. Vertical gradients in velocity within the ice column are neglected when plotting Fig. 3 and fig. S3. This is justified given the low deformation rates recorded at the field location.

Modeling experiments

A suite of flowline models for a slab glacier on an inclined plane were implemented using Elmer/Ice (15). The model solves the Navier-Stokes equations in two dimensions for an incompressible fluid with a nonlinear, temperature-dependent, Glen’s style rheology. The slab glacier has a thickness of 650 m, has a length of 800 m, and is resting on a bed with a mean slope of 1°. The model was run on a mesh with an approximate resolution of ~1.5 m generated by ElmerGrid. The lateral flow boundaries of the model are periodic. The surface slope is initially set to 1° and was allowed to evolve as a free surface until equilibrium was reached.

The model assimilates ice temperatures collected at the center of the borehole grid and thus includes a realistic representation of vertical rheologic variations. To control the sliding rate, a linear Weertman-style sliding law was prescribed at the basal boundary that relates the local basal shear stress to the sliding velocityub=τbC(6)where ub is the basal sliding velocity, C is the slip coefficient, and τb is the local basal shear stress. A Weertman-style sliding law is a simplification of the complex physics that operate at the bed. However, the goal of the modeling is to elucidate the near-basal ice dynamical response of sliding over small-scale roughness, making it sufficient for prescribing sliding velocities at the basal boundary.

Different modeling scenarios were tested in an effort to reproduce the main tilt rate characteristics in the dataset. In each model scenario, ice was forced to slide over simplified, yet plausible, basal boundaries. A full list of model scenarios tested is cataloged in table S2 and split into two broad categories: roughness and patchiness. Sliding over basal roughness was modeled for three different types of basal topography: sinusoidal roughness, asymmetric roughness, and bed asperities. The sinusoidal roughness and asymmetric roughness test how ice flows over smooth repeating bedrock roughness with a simplified but plausible shape. Bed asperities were modeled as small-amplitude Gaussian bumps with spacing equal to five times the bump width across the bed.

Basal patchiness was modeled by imposing alternating sticky and slippery patches along a smooth basal boundary. The slip coefficient of slippery patches was set to a value about two orders of magnitude more slippery than the sticky patches to mimic a strong contrast in slipperiness. We also considered a scenario that combined roughness and patchiness. Motivated by the common conceptualization of linked cavity basal hydrology, we coupled the sinusoidal bed and patchy bed scenario (termed “sine slip”) by imposing slippery patches on the lee side of bedrock features.

A total of 13 models were run each with a unique basal boundary. We limit the roughness and patchiness to wavelengths of 80 to 100 m because we are interested in replicating the characteristics of longest wavelength tilt rate variations observed. The sliding coefficient for each model run was tuned until the surface velocities approximated those measured at the field location (114 ± 10 m a−1) averaged over the 2015–2016 winter. A table of model parameters for all model scenarios is found in table S2.

Comparison of model results to basal tilt rates

For each model scenario, synthetic tilt rates along theoretical inclinometer flow paths were compared to the observations (Fig. 3 and fig. S5). The theoretical flow paths were calculated using the modeled velocity field and the installation height above the bed as the starting location. In general, the flow path of the inclinometer mimics the basal topography below it, but with decreasing amplitude moving up through the ice column. Roughness of the bedrock topography means that different flow paths exist depending on whether the inclinometer is installed above a peak and trough of the bed roughness. However, we found very little difference between the theoretical inclinometer flow paths starting from either the trough or ridge given the same installation height.

Synthetic tilt rates were calculated along the flow paths from the gradients in the velocity field, yielding a synthetic record of inclinometer tilt rates versus basal distance to compare to the observations. It is conceivable that inclinometers could be located above or below their installation depth due to uncertainty associated with cable stretch or ice flow before the 2015–2016 winter. This uncertainty was taken into account when finding the boundary best fit to the observations by considering synthetic tilt rate records from flow paths starting at ±2 m the installation height. An example of modeled inclinometer tilt rates for each type of basal boundary is plotted in fig. S5. The observations were compared to both the synthetic tilt rate records and the spatial pattern of tilting. The best fit to the observations is plotted in Fig. 3.

Sensitivity of deformation rates to rheologic uncertainty

To investigate the range of possible deformation rates that could be estimated if the rheology was unconstrained by in situ measurements, we used the driving stress and a range of three common rheologic parameters to estimate the deformation velocity for the field location (fig. S8).

The parameters tested were as follows: Enhancement factor: Pleistocene age ice softens and enhances the rate of deformation (24). Enhancement factors, E, have been estimated at E ≈ 2 (9), E ≈ 4 (10), and E ≈ 10 (16, 17) for Pleistocene ice in Greenland. Near the margins in Greenland, Pleistocene ice is found to make up the bottom 15 to 20% of the ice column (39). We applied the above enhancement factors to the bottom 20% of the ice column. Temperature: Ice temperature is the primary control on the creep parameter A. Warmer ice deforms more quickly than colder ice (24), and thick temperate layers can form near the margin (27, 28) where water content can soften the ice further (24). We calculated the creep parameter using uniformly warmed (+3° and +5°C) temperature profiles from that measured at the field site. Flow exponent: While n = 3 is the most commonly used flow exponent, the deformation velocity was calculated using n = 4, which has also been suggested for Greenland ice (14, 40). Other rheologic factors such as ice fabric, grain size, and water content also influence rheology (24). These were not specifically quantified because of limited constraints on rate enhancement of such factors and lack of knowledge on their occurrence and distribution on the GrIS. Even so, the influence of some of these factors is already considered in the enhancement factors tested.

Computation of driving stress across the Western GrIS’ ablation zone

Driving stress across western GrIS’ ablation zone was calculated from Eq. 3 using a combination of CryoSAT ice sheet surface elevation (41) and bedrock topography (42). Ice thickness was computed by differencing the surface elevation and bed elevation datasets. Surface slopes were computed by smoothing the surface elevation using a weighted average incorporating neighboring grid cells. Gradients of the smoothed product were then computed using a centered difference with 2-km spacing.

Supplementary Materials

Supplementary material for this article is available at

Fig. S1. Radiogram of field location.

Fig. S2. Deformation profile summary.

Fig. S3. Basal tilt rate versus basal distance.

Fig. S4. Basal tilt rate, surface velocity, and temperature.

Fig. S5. Modeled basal tilt rates.

Fig. S6. Modeled deformation profiles.

Fig. S7. Near-basal stress field.

Fig. S8. Deformation velocity profiles for different rheologic parameters.

Table S1. Inclinometry data summary.

Table S2. Flowline model parameters.

Supplementary analysis

Reference (44)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This project was funded by the NSF Office of Polar Programs-Arctic Natural Sciences awards #1203451 and #0909495, SKB, NWMO, Posiva Oy, and NAGRA. Author contributions: N.M. analyzed the data, produced Figs. 1 to 3, and is the primary author of the manuscript. N.H., J.H., and T.M. all contributed significantly to the writing and editing of the manuscript. T.M. produced Figs. 4 and 5. N.H. and J.H. designed the field experiment. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. All data used to generate this manuscript will be archived at

Stay Connected to Science Advances

Navigate This Article