Login using

You can login by using one of your existing accounts.

We will be provided with an authorization token (please note: passwords are not shared with us) and will sync your accounts for you. This means that you will not need to remember your user name and password in the future and you will be able to login with the account you choose to sync, with the click of a button.

1Atmospheric and Oceanic Sciences, University of Maryland, College Park, College Park, MD, United States

2Geoscience and Mathematics, Penn State DuBois, DuBois, PA, United States

3Jet Propulsion Laboratory, La Cañada Flintridge, CA, United States

The Greenland Ice Sheet has experienced accelerated mass loss over the last couple decades, in part due to destabilization of marine-terminating outlet glaciers. Retreat and acceleration of outlet glaciers coincides with atmospheric and oceanic warming resulting in a significant contribution to sea-level rise. The relative role of surface meltwater production, runoff and infiltration on the dynamics of these systems is not well-understood. To assess how surface meltwater impacts shear margin dynamics and regional ice flow of outlet glaciers, we investigate the impact of basal lubrication of Jakobshavn Isbræ shear margins due to drainage from water-filled crevasses. We map the areal extent of inundated crevasses during summer (May–August) from 2000 to 2012 using satellite imagery and determined an increasing trend in the total areal extent over this time interval. We use a numerical ice flow model to quantify the potential impact of weakened shear margins due to surface melt derived basal lubrication on regional flow velocities. Ice flow velocities 10 km from the lateral margins of Jakobshavn were amplified by as much as 20%, resulting in an increase of ~0.6 Gt yr−1 in ice-mass discharge through the shear margins into the ice stream. Under future warming scenarios with increased surface melt ponding, simulations indicate up to a 30% increase in extra-marginal ice flow. We conclude that surface meltwater will likely play an important role in the evolving dynamics of glacier shear margins and the future mass flux through Greenland's major marine-terminating outlet glaciers.

Introduction and Background

Over the last decade, the Greenland Ice Sheet has experienced accelerated mass loss resulting in a doubling of its contribution to sea-level rise (Shepherd et al., 2012). Previous work asserts that the majority of Greenland's mass loss is due to solid ice discharge from marine-terminating outlet glaciers (Zwally et al., 2011; Joughin et al., 2012, 2014; Moon et al., 2012). Much of the mass loss is associated with observed increases in velocity in the lower trunks of Greenland's major outlet glaciers, which has been attributed to oceanic forcing of glacier termini through incursion of warm water into fjords (Joughin et al., 2012; Walter et al., 2012; Bondzio et al., 2016; Muresan et al., 2016). Recently, the partitioning of mass loss has been reassessed such that ~60% is due to changes in surface mass balance and the remainder (~40%) from solid ice discharge since 1991 (van den Broeke et al., 2016). Infiltrated surface meltwater has generally been implicated in enhancing ice-shelf calving instability and lubrication of previously frozen ground (Parizek and Alley, 2004; Shannon et al., 2013; Tedstone et al., 2013; Everett et al., 2016). Although surface melt may be a less important control on dynamic acceleration over Greenland's outlet glaciers than changes in glacier force balance initiated at the ice-ocean interface (Vieli and Nick, 2011), the effects of meltwater change on ice flow into the fast-flowing outlet glaciers has yet to be fully investigated.

In this analysis, we focus on Jakobshavn Isbræ, one of Greenland's major fast-flowing outlet glaciers, which like many is confined to a deeply incised trough that promotes well-established shear margins. Despite its narrow trough, high driving stresses on Jakobshavn impel large and, over the last three decades, accelerating ice fluxes that are currently draining ~7% of the ice sheet (Bindschadler, 1984; Echelmeyer et al., 1991; Rignot and Kanagaratnam, 2006; Cuffey and Paterson, 2011; Moon et al., 2012). Near Jakobshavn's terminus, lateral drag compensates ~40% of the gravitational driving stress responsible for mass discharge, while upstream basal drag accommodates ~20% (Van der Veen et al., 2011).

Reduction in lateral drag has been implicated in observed changes in Jakobshavn ice flow over the last couple of decades (Van der Veen et al., 2011). Vieli and Nick (2011) suggest acceleration of Jakobshavn was sustained from 1998 to 2001 due to enhanced mechanical weakening of the shear margins via strain heating despite a slow-down in terminal retreat. Earlier work by Lampkin et al. (2013) examined drainage of water-filled crevasses (called “saturated crevasses”; CV) within the shear margins of Jakobshavn Isbræ during the 2007 melt season and estimated meltwater drainage volumes (~9.23 × 10−3 ± 2.15 × 10−8 km3) on the order of the large supraglacial lakes found outside of the shear margins. They proposed that drainage of water-filled surface crevasses could lead to enhanced meltwater-induced basal sliding beneath the shear margins, thereby promoting a so-called “hydrologic weakening” of the shear margins that could accelerate ice flow and mass discharge both within and outside of the ice stream.

Furthermore, Cavanagh et al. (2017) documented an increase in extra-marginal ice flow of ~20% associated with observed drainage from water-filled crevasse ponds during the 2013–2015 melt seasons and consistent with numerical simulations that predict a potential increase in flow of ~35%. Joseph and Lampkin (2017) examined the temporal variability of hydrologic state (drained or filled) for each crevasse system using a fusion of data from several optical sensors and found significant annual variability in duration and timing of drainage. They also observed that some systems demonstrated multiple drainage events within a single melt season. Recently, Bondzio et al. (2016) utilized a thermomechanical ice flow model to assert that calving-induced weakening of the shear margins (through changes in the englacial stress state that lead to a reduction in effective viscosity) can account for widespread increases in regional ice flow between 5 and 10%. Their work supports the earlier assertion by Vieli and Nick (2011) that rheological damage may be an important process contributing to Jakobshavn's acceleration.

Thus, several factors likely play important roles in influencing contemporary and future changes in shear-margin dynamics. Each of these processes is likely influential over a range of spatial and temporal scales and varies as a function of changes in external forcing. The impact of each process on regional ice flow and their interactions are not well-understood.

Objectives

Given insights from previous efforts, it is likely that several processes will impact the shear margins in a changing climate. We identify three distinct processes: (1) strain softening, (2) basal lubrication from pond drainage and infiltration (hydrologic shear weakening), and (3) cryo-hydrologic warming (Figure 1). Strain softening results from the dissipation of heat due to deformation of the ice (Cuffey and Paterson, 2011) driven by lateral drag and vertical shearing in the shear margins (Harrison et al., 1998; Whillans and van der Veen, 2001; Van der Veen et al., 2011; Lüthi et al., 2015; Poinar, 2015). Heat production via lateral drag can be large and is proportional to the strain rate and lateral shear stress (Whillans and van der Veen, 2001). Cryo-hydrologic warming of the ice column can occur due to freezing of englacial melt water in deep crevasses and racture networks, resulting in heating due to latent heat flux (Phillips et al., 2010, 2013; Colgan et al., 2016; Poinar et al., 2017) such that ice temperatures can potentially increase by 1.8°C for 1% water retention per ice volume (Colgan et al., 2016). The relative impact of each of these processes on regional ice flow, their spatial and temporal variability and feedbacks between them is not well-understood. Because of the importance of lateral drag in Jakobshavn's force balance, we will focus solely on hydrologic shear weakening through a first-order investigation into the impact of basal lubrication of Jakobshavn Isbræ shear margins due to meltwater drainage from water-filled crevasses. We focus on this process because the emergence of ponded water within the shear margins is a significant development for Greenland's outlet glaciers. Additionally, the delineation of these ponds to assess changes in their properties and the implementation of an end-member parameterization of their influence on ice flow are easily implemented using visible satellite imagery and numerical ice-flow models. We map the areal extent of inundated crevasses during summer (May–August) from 2000 to 2012 using satellite imagery and determined an increasing trend in the total areal extent over this time interval. We use a numerical ice flow model to quantify the potential impact of weakened shear margins due to surface melt derived basal lubrication on regional flow velocities. In this effort, modeling scenarios provide an estimate on the potential impact of meltwater on ice flow due to hydrologically weakened shear margins under a warmer regional climate. In this analysis, we evaluate the hypothesis that weakened shear margins could accelerate ice discharge into the ice stream. Therefore, we focus on how extra-marginal flow is impacted. For a full summary of symbols and variables used in this analysis see Table A1.

FIGURE 1

Figure 1. Conceptual graphic depicting major processes within the shear margins that would become activated in a warming climate. Processes include strain heating enhanced by acceleration of ice (calving induced), basal lubrication, and cryo-hydrologic warming due to infiltration of surface melt from water-filled crevasse drainage and surface melt runoff.

Data and Methods

Water-Filled Crevasse Mapping

The imagery used for this analysis are Landsat-7, Enhanced Thematic Mapper Plus (ETM+) Scan-Line Corrector (SLC)-off data acquired from the US Geological Survey LPDAAC. Landsat-7 panchromatic (Band 8) images have a spatial resolution of 15 m. Cloud-free, Landsat scenes were acquired over the May, June, July, August, and September months from 2000 to 2013 (https://landsat.usgs.gov/landsat-data-access). Many studies have utilized several techniques to map supraglacial features (Box and Ski, 2007; McMillan et al., 2007; Sneed and Hamilton, 2007; Sundal et al., 2009; Selmes et al., 2011; Banwell et al., 2014; Leeson et al., 2015; Pope et al., 2016; Williamson et al., 2017; and several others) using both manual and automated routines. Most automated routines rely on spectral information to delineate ice from water to determine the extent and depth of water in supraglacial lakes. Automated and manual techniques have many advantages but are both limited by the spatial resolution of specific satellite systems, cloud cover, heterogeneity in surface and illumination conditions (Sneed and Hamilton, 2007; Leeson et al., 2015; Williamson et al., 2017). The water-filled crevasse systems examined in this analysis occur in rough terrain and are susceptible to reduced reflectance due to shadowing. Because this degree of shadowing limits the utility of automated routines, we chose to manually retrieve the areal extent of these features. The boundaries of each water-filled crevasse region were manually digitized by a trained analyst based on the observable difference in reflectance between water and ice. A second trained analyst was used to review mapped boundaries as quality control. Errors in mapped boundaries are a result of the spatial resolution of the imagery and instances where missing scan lines cut across features. Missing scan-lines repeat at regular intervals throughout the scene with increasing width near the scene edges. Scene centers contain regions with few missing scan-lines and were usually 1–2 pixels in width, whereas those near the edge were about 15 pixels in width on average (Lampkin, 2011). Most water-filled crevasse features occupied regions near the center of the scene. Scenes where water-filled crevasse ponds were near the edge, where missing scan lines substantially obscure the ability to visually interpolate pond boundaries were not mapped. Scenes with fractional cloud cover that substantially obscured a feature were also omitted. An estimate of absolute uncertainty in delineated area due to the limits of spatial resolving capability of Landsat (γlandsat)and missing scan-lines (γscanline) for the best case (2 scan-line pixels) given by (γlandsat)2+(γscanline)2 is 2.25 × 10−4 km2 and for the worst case (15 scan-line pixels) is 3.18 × 10−4 km2. Absolute error for any given mapped feature could vary between these two values but the overall error is closer to the best case value given most water-filled crevasse groups were mapped near the center of a scene. We solely utilized Landsat imagery to map the areal extent of each water-filled crevasse group from 2000 to 2012 (Figure 2A) and do not leverage other visible imagery. We are primarily interested in characterizing the trend in areal extent at the annual scale and not seasonal variability. The Landsat archive has sufficient samples to determine the trend in extent over the study period. Joseph and Lampkin (2017) have evaluated the seasonal variability in hydrologic state (filled vs. drained) and additional work will be completed to expand this assessment to changes in inundation extent.

Surface Temperature

Data from the Program for Regional Climate Assessment (PARCA) GC-Net archive (Steffen et al., 1996) were used to examine the temporal variability in near surface air temperatures (http://cires1.colorado.edu/steffen/gcnet/). Monthly mean, minimum, maximum, and standard deviations of air temperatures between June and August of each season were derived from a combination of data from the JAR1 and Swiss Camp meteorological stations (which are ~38 and 50 km north of Jakobshavn, respectively) from 2001 to 2013 (Table 1).

TABLE 1

Table 1. GC-Net station locations.

Modeling of Extra-Marginal Ice Acceleration Due to Hydrologic Shear Weakening

Model Description

Basal lubrication of Jakobshavn shear margins is expected to have two consequences: speed-up within the main trunk due to reduced lateral shear, and acceleration of ice outside of the main ice stream resulting in expansion of the catchment. In this analysis, we investigate the extent to which drainage of water-filled crevasses could enhance basal lubrication and, in turn, the impact of this basal lubrication on ice velocities and discharge from regions beyond the shear margins of the outlet glacier. Surface meltwater accumulates among assemblages of crevasses within the shear margins (Figures 2A,B). Drainage of water-filled crevasses occurs and can result in lubrication of the bedrock/ice interface. Here, we assume that drained meltwater infiltrates completely to the bed. Enhanced mass flux through the margins is a result of a melt-induced decrease in shear arising from increased basal sliding. This mass flux perturbation is transmitted upstream due to longitudinal stretching, increasing along-flow ice motion and upstream mass flux (Lampkin et al., 2013) (Figure 2C). This process is similar to that described by Christoffersen et al. (2018), where lake drainages can induce other drainage events tens of kilometers away. Reduction in basal traction due to infiltrated melt water enhances membrane stresses resulting in regional ice acceleration and tensile shock (Christoffersen et al., 2018). We hypothesize that a similar process can result from melt water injected into the shear margins of Jakobshavn. This impact is most extensive along the southern lateral margin of Jakobshavn and therefore we focus our analysis on this region.

For clarity, we do not explicitly consider the influence of cryo-hydrologic warming or strain softening on lateral shearing within the margins. Refrozen meltwater within the ice column is a viable mechanism to weaken the shear margins through cryo-hydrologic warming of the ice column as latent heat within infiltrated water is released during refreezing (Van der Veen et al., 2011). Cryo-hydrologic warming has not been quantified in Jakobshavn but has been assessed in the ablation zone of Sermeq Avannarleq (~75 km north of Jakobshavn), where the upper 100 m is estimated to warm by ~10 K for an average spacing between englacial fractures (R) of ~20 m (Phillips et al., 2010). The rate of warming in the ice column increases when R is small (Phillips et al., 2010). Given that crevasses along the shear margins of Jakobshavn are densely spaced, we might expect that infiltrated meltwater may be effective in reducing the ice viscosity through warming in locations where water can readily enter without drainage to the bed. More work is necessary to thoroughly investigate this process. Given we are specifically concerned with evaluating the hypothesis that hydrologic shear weakening can amplify extra-marginal ice flow (Lampkin et al., 2013), we do not comprehensively investigate the feedbacks between basal lubrication of the shear margin and ice flow within the main ice stream, but rather provide a cursory examination of this component of the outlet glacier system to assess whether it could influence its flow and future evolution.

First-order assessment of the impact of drainage on regional ice flow was accomplished through the use of ISSM (Larour et al., 2012) for the development of 2-D plan-view diagnostic and prognostic models. Ice is treated as an isotropic, viscous material, which includes the incompressibility and momentum balance equations. Assuming negligible acceleration and neglecting the Coriolis Effect, the equations for conservation of mass and momentum are, respectively:

∇·v=0(1)

∇·σ+ρig=0(2)

where v is velocity, ρi ice density, g gravitational acceleration, and ∇ • σ is the divergence of the stress tensor σ. Using a nonlinear constitutive flow law for ice (Glen, 1952), the effective ice viscosity, μ, is defined as:

2μ=Bε.e1-nn(3)

where n is 3, B is the ice hardness parameter, and εe. is the effective strain rate.

We use the shelfy-stream approximation to the Stokes equation to model ice flow (Morland, 1987; MacAyeal, 1989), which simplifies the full-Stokes stress balance to a system of two equations with two unknowns, under the assumption that vertical shear and bridging effects are negligible.

The model domain encapsulates all 7 water-filled crevasse groups (CV1-7) (Figure 2A) and was populated with an unstructured triangular mesh with uniform resolution set at 200 m. Surface and basal topography were prescribed using a 150-m resolution BedMachine DEM with a nominal date of 2007 (Morlighem et al., 2014). Gridded geothermal heat flux and surface mass balance values used in this analysis were derived from the Sea-level Response to Ice Sheet Evolution (SeaRISE) common model data archive (e.g., Shapiro and Ritzwoller, 2004; Ettema et al., 2009; Bindschadler et al., 2013; Nowicki et al., 2013). Ice thickness is calculated as the difference between surface and basal topography. Surface velocity was prescribed based on TerraSAR-X data from the MEaSUREs Greenland Ice Sheet velocity, select glaciers archive (Joughin et al., 2010a). Surface velocity component fields (u, v) at a nominal spatial resolution of 100 m were downloaded from the NSIDC MEaSUREs data archive for the Jakobshavn focal region (http://nsidc.org/data/nsidc-0481/). The model only used velocity data from the summer of 2009 as a representative melt season. The 2009 melt season provided the best spatial coverage of the extra-marginal ice field in the vicinity of the water-filled crevasse ponds in addition to being one of several seasons with anomalous speeds ups before 2009 (Joughin et al., 2014). Models were derived without basal hydrology and thermal solution schemes, resulting in depth-integrated temperatures remaining constant. The upper boundary of the model is considered stress-free. A modified Coulomb-style basal friction (τb) law:

τb=-k2Nr||vb||s-1vb(4)

is used at the ice-bedrock interface (Van der Veen et al., 2011), where vb is basal velocity, N effective pressure, s and r are friction law exponents, and k is the friction coefficient. Friction accounts for the impact of a thin layer of water at the bedrock/ice interface via N:

N={g(ρiH+ρwb),b≤0,w=0ρigH,b>0,w=0fiρgH,∀b,w>0(5)

This is defined for locations where the bedrock elevation is above sea-level (b > 0) or below (b < 0) as a function of ice density, water density (ρw), water layer thickness at the base (w), and ice thickness (H). Where a thin water layer is present, a pre-factor (f) is included to account for the proportion of glaciostatic pressure that is compensated by hydrostatic pressure. The modified basal friction scheme in the model scenarios is configured to be invariant to the magnitude of basal water layer thickness, but is sensitive to the degree to which hydrostatic pressure compensates the glaciostatic pressure.

Model Scenarios

Simulations were executed in three phases; basal friction inversion, diagnostic model simulating velocity field in the absence of a water layer at the bed (baseline or reference model) and a suite of diagnostic models with varying degrees of basal lubrication at the bed due to drainage from water-filled crevasses (perturbation models). Phase I establishes a reference basal friction field from which lubrication enhancements due to infiltrated meltwater are applied. This reference friction field represents the current state in basal friction based on a quantitative inversion of the basal friction coefficient field. Numerical inversion is implemented using a control method (MacAyeal, 1983; Larour et al., 2005), which finds the distribution in basal friction through gradient minimization of a cost function (J), which is a measure of the misfit between observed (vx, vy) and modeled (υxm, υym) surface velocity components. The cost function is given by:

J(υ)=∬S12((υym−υy)2+(υxm−υx)2)dS(6)

where the criteria for simulation termination is when J reaches a minimum and variability is < 1%. The estimated k field was then used across the bedrock/ice interface for the following diagnostic models.

Two sets of model simulations were implemented to evaluate the relative impact of shear margin lubrication: (i) a diagnostic baseline model (Phase II) without lubrication that provides estimates of velocity (vphaseII) when the w-field within the composite extent of water-filled crevasse regions is set to (0) and the inverted basal friction coefficient (k) from Phase I is used throughout the domain and (ii) another set of diagnostic and prognostic model runs with varying degrees of basal lubrication that (Phase III).

Phase III simulations provide an estimate of velocity (vphaseIII) resulting from complete vertical infiltration of surface water and basal lubrication within the composite extent of water-filled crevasses (w > 0). These scenarios use the modified Coulomb style basal friction scheme (Equations 4 and 5). An indication of the impact of drainage from water-filled crevasses on regional ice flow is quantified through the relative difference metric given by:

Δv=(vphaseII-vphaseIII)/vcv¯(7)

which normalizes the difference between Phase II and III modeled velocities expressed as percentages, where (vcv¯) is the mean between modeled Phase II and III velocity estimates. This metric is used to evaluate both diagnostic and prognostic output arising from the perturbations detailed below.

Description of Experiments

Diagnostic models explore the impact of water-filled crevasses drainage on extra-marginal ice flow through varying basal hydraulic pressure (Δνdf) and areal extent (ΔυdArea) of crevasse ponds using Equation (7). Changes in these parameters are evaluated along a flowline through the southern margin of the ice stream. The initial value for f was set such that the basal hydraulic pressure compensates for 99% of the glaciostatic pressure beneath all water-filled crevasse groups based on bore-hole measurements (Lüthi et al., 2002). Such a high basal hydraulic pressure would be indicative of an inefficient subglacial hydrologic network within a distributed configuration or a coupled linked-cavity conduit system with under-developed channels forced by short duration meltwater input (Schoof, 2010; Flowers, 2015). Within this set of Phase III experiments, we examine ice flow sensitivity to f-values of 25, 50, and 75%. Diagnostic model sensitivity to varying areal extent by 50% of the decadal composite extent (f = 99.9%) was implemented as well.

Depth-average mass flux due to hydrologic shear weakening through a fluxgate along the southern margin (Figure 2A) was determined for Phase II and III (for f = 99.9%) diagnostic models. The southern margin fluxgate is ~23 km in length and is parallel to the dominant direction of ice flow and is located ~3 km away from the centerline. Total mass flux through the gate is determined by:

Mi=∫s=0s=LρiceHυ¯·ndl(8)

where s is the position along the curvilinear fluxgate boundary, L horizontal length of fluxgate, ρice depth-average ice density over the ice thickness along the fluxgate, υ¯ depth-average velocity vector, n local normal to fluxgate contour, and dl differential curvilinear distance along the fluxgate.

Phase III prognostic simulations were developed to quantify ice flow response to hydrologic shear margin weakening due to seasonal forcing in basal pressurization, areal extent of lubrication, and lubrication duration. Transient solutions for these models were run at the same mesh resolution as the diagnostic models. A mass conservation solution was included, which estimates the rate of ice thickness change as a function of surface mass balance (Ṁs), basal mass balance (Ṁb), and ice flow divergence given by:

∂H∂t=Ṁs+Ṁb-∇·Hv¯(9)

where H is ice thickness, v¯ depth-average horizontal velocity. The shelfy-stream ice flow model was used as described above. Transient solutions were run for 10 years from 2017 to 2027 using surface mass balance forcing from 2009. We used a constant surface mass balance to isolate the impact of basal lubrication. A grid-independence analysis was conducted to determine a stable timestep size for a fixed mesh resolution which satisfies the Courant-Friedrichs-Lewy stability criterion (Courant et al., 1928). A 2-week time step was determined to meet the stability criterion. Transient baseline and perturbation models were forced with the same boundary, velocity, and surface forcing data described above. Because we are interested in the relative difference between the various basal lubrication models and the baseline (no lubrication) scenario, model spin-up is initialized with the diagnostic model fields under the assumption of no additional basal lubrication.

Normalized differences in extra-marginal velocity between transient baseline and perturbation scenario is determined for varying basal pressurization (υtf), areal extent (ΔυtArea), and lubrication duration (Δυtduration) using Equation (7). At each time step, the maximum value of these perturbation metrics along the analysis flowline are derived to spatially summarize the impact of hydrologic shear weakening. υtf is computed for f parameters of 99.9, 50, and 10%. Estimates of ΔυtArea are quantified, where the composite areal extent across all CV groups is uniformly scaled by annual deviations from a reference temperature (ΔTA). ΔTA is set to 4.7°C based on averaging 2-meter July air temperatures from 2004 to 2016 (Ettema et al., 2009) over the water-filled crevasse regions. The month of July is used to determine ΔTA because this month corresponds to a maximum in both summer-time monthly mean air temperatures and the peak areal extent of ponding (Lampkin et al., 2013).

During transient runs, the duration of basal lubrication is extended beyond the 30 days used in steady-state experiments to 60 and 90 days. These intervals are in 30-day increments to bracket the potential time it takes for basal meltwater to evacuate the local region where surface melt drained to the bed. The rate of subglacial drainage has not been established throughout the study region, but the drainage rate and subglacial topography among other variables are important (Werder et al., 2013). Short duration drainage from the water-filled crevasse systems would likely support a distributed subglacial hydrologic network, which is less efficient than a conduit (Schoof, 2010; Werder et al., 2013). Cavanagh et al. (2017) examined subglacial hydraulic potential throughout Jakobshavn's shear margins in the vicinity of the water-filled crevasses and speculated that some crevasse ponds may develop efficient channelization due to relatively large hydraulic potential gradients. This would result in rapid evacuation of subglacial water and lower lubrication durations. Because the focus of this study is on the general impact of basal lubrication within the shear margins, we do not provide a comprehensive treatment of the complex subglacial hydrologic system. The implications of this limitation are discussed later.

Results

In this work, water-filled crevasses were mapped as spatially coherent groups (CV groups in Figure 2A). The composite areal extent of inundated regions within the shear margins over the study interval comprises a substantial fraction of the shear margins but is spatially limited to locations with small inflections in ice-flow direction (Figure 2A, blue regions). The ponded areas are spatially consistent with basal topographic basins, which likely control where they occur each season (Lampkin et al., 2013). The total area summed from 2000 to 2012 shows CV 6 and 7 to be the largest, while CV 1 through 5 are less extensive (Table 2). Over the study period, mean summer (June–August) surface temperatures have increased from ~1°C in 2001 to ~2.5°C in 2012, with cooler mean summer temperatures slightly below freezing in 2011 and 2013 (Figure 3). Although we are not able to capture inter-seasonal and seasonal variability in areal extent of the water-filled crevasse groups due to cloud cover in the imagery, the number of mapped scenes increases along with mean monthly surface temperatures, likely resulting in capturing near-maximum areal extents (Figure 4A). We derive seasonal total ponding area (sum of area of all water-filled groups for a given summer) and a weighted seasonal total area based on scene acquisition date and their departures from mean peak summer temperatures. Total areal extents demonstrate an increasing trend over the study period at ~0.32 and ~0.35km2/year for unweighted and weighted cases, respectively (Figure 4B).

TABLE 2

Table 2. Total composited area of water-filled crevasses from 2000 to 2012 using Landsat visible imagery.

FIGURE 3

Figure 3. Time series of near surface air temperatures during summer months (June–August) from a data set combined from GC-Net meteorological stations JAR1 and Swiss Camp from 2001 to 2013. Time series show summertime mean, maximum, minimum, and standard deviation of air temperatures over the study period.

FIGURE 4

Figure 4. Histograms of mapped image dates and monthly air temperatures. (A) The number of image dates used to map water-filled crevasses per month (left y-axis). Mean monthly near surface temperatures measured from JAR1 and Swiss Camp stations (right y-axis). (B) Total areal extent of water-filled crevasse groups from May to September for each year. Areal extent trends shows an increase over study period for unweighted (solid line) extent. A weighted total area was computed (dotted line) by fitting a Gaussian to the image histogram in (A) and calculating annual total areal extent weighted by the month any particular water-filled boundary was mapped.

We quantify the impact of basal lubrication of the shear margins on ice flow as the relative difference in modeled velocities between the baseline water-free model and those with basal lubrication (Figure 5A). Velocity differences show variability both within and outside of the main trough of Jakobshavn (Figure 5A). The largest increase in velocity is concentrated in the extra-marginal ice field near the zones of basal lubrication at significant distance beyond the shear margins. Faster flow due to basal lubrication of the margins is spatially more extensive in the southern extra-marginal field, while regions of similar speed-up are more concentrated parallel to the northern margin within a distance of ~6 km (Figure 5A).

FIGURE 5

Figure 5. Impact of hydrologic scenarios for shear weakening on diagnostic ice velocity. (A) Basin-wide percentage change in modeled velocity for Phase II and III diagnostic scenarios. Drainage results in an increase of up to 20% in extra-margin ice speeds extending as far as 10 km from the shear margins and an increase of about 5–7% within the trough. (B) Percent change between diagnostic Phase II and III velocities along our flowline (red line in Figure 1A) for varying values of f (proportion of glaciostatic pressure compensated by basal water pressure) and a 50% increase in areal extent of all CV groups. Arrows indicate ice-flow direction.

Sensitivity tests were designed to evaluate the impact of varying the magnitude of basal pressurization (defined by the ratio of hydrostatic to glaciostatic pressure, f) and the lubrication area (A). Along our flowline (Figure 2A), as expected, percentage increase in diagnostic velocities is linear with pressurization (Figure 5B). The maximum ice accelerations of ~10% (f = 99.9%) and ~1% (f = 25%) occur ~3 km from the main trough. Furthermore, uniformly increasing CV area by 50% while holding pressurization at 99.9% leads to an additional 6% increase in ice speed (on par with the simulated impact of f = 75%), which corresponds to a mass flux through the southern margin of ~16 Gt a−1 and also results in enhanced flow over a much larger region (Figure 5B). The areal extent of drainage thereby has a greater influence on regional ice flux than the magnitude of basal pressurization.

We spatially evaluate the impact of water-filled crevasse drainage on extra-marginal ice flow through derivation of zonal spatial means in velocity differences between the water (f = 99.9%) and no-water diagnostic model cases (Δv¯zone). The zones buffer the center of the ice stream at 1 km increments (Figure 6A). Δv¯zoneis ~ 6% near the center of the glacier, reaches a maximum of 20% at about 4 km beyond the shear margins and decreases to about 1% at 10 km. The range in Δv¯zonedecreases commensurately with distance from the margin (Figure 6B). Within the trough, diagnostic results indicate that drainage from water-filled crevasses can increase ice flow on average by ~5–7%, with a maximum of 12%.

FIGURE 6

Figure 6. (A) Map depicting ten nested zones (gray lines) at 1 km increments extending 10 km away from center of the outlet glacier (green line) within the study domain (red box). These zones were used to quantify spatial variability of Diagnostic model results for the impact of hydrologic shear weakening due to drainage from water-filled crevasses (Main text; Figure 3A). Superimposed are the decadal composited extent of water-filled crevasse systems (blue) used in this analysis. (B) Mean, maximum, minimum, and standard deviations of percentage change in velocity between the diagnostic baseline model and basal lubrication scenario for f = 99.9% within each zone in (A).

We further explore the impact of hydrologic shear weakening on extra-marginal flow from 2017 to 2027, during which we explore sensitivity to lubrication duration, inundation area, and basal pressurization. Transient runs indicate that lubrication durations at 60 and 90 days result in slightly larger speed up (~16%) (Figure 7A) than 30 days (Figure 7C). We vary the inundation area as a function of summer-time (July) temperature deviations from ΔTA. A maximum increase in extra-marginal ice speeds of ~30% is simulated when the areal extent is increased by ~13% (Figure 7B). Additionally, the speed-up is larger as f increases for a transient forcing of 30 days, which is similar to our results from the Diagnostic model (Figure 7C). Interestingly, the time for extra-marginal speeds to return to pre-perturbation magnitudes is longer for smaller f. This is likely the result of small, post-lubrication driving stresses (vis-à-vis smaller surface slopes), which slowly flush extra-marginal ice back into the margins.

FIGURE 7

Figure 7. Phase III transient model scenarios of ice velocity due to hydrologic shear weakening from 2017 to 2027. (A) Maximum percentage change in velocity (left y-axis) along our flowline (Figure 1A) between the transient baseline model and a temporally varying lubrication duration at 60 and 90 days for f set to 99.9%. (B) Maximum percentage change in velocity along our flowline between the transient baseline model and a temporally varying extent in ponded area (left y-axis). Magnitude of areal extent dilation proportional to the deviation of annual July temperatures above the reference temperature (4.7°C) (green points, right y-axis). (C) Maximum percentage change in velocity along our flowline between the transient baseline model and scenarios for varying values of f.

Discussion

In this analysis, we assume hydrologic shear weakening is due solely to basal sliding. However, it is unlikely that all CV groups weaken the shear margins to the same degree. The larger, spatially more cohesive CV groups, which are found in the lower reaches of the ice stream where ice is thinner, are more likely to deliver water directly to the bed. In contrast, the CV groups (6 and 7) that are defined by isolated ponded regions that are distributed over a larger area at higher elevations where the ice is thicker, are less likely to instantaneously drain large volumes of water directly to the bed (Figure 2A). The capacity to propagate water-filled fractures to the bed is dependent on the rate of meltwater input into the fractures, which controls the degree to which hydrostatic pressure can overcome fracture tip closure (van der Veen, 1998; Van der Veen, 2007; Alley et al., 2005). Though the water-filled crevasses in this study are closely spaced, there is sufficient tensile stress over these features to propagate fractures to the bed (Lampkin et al., 2013). Given that water-filled crevasse groups largely occur within local depressions (Lampkin et al., 2013) where there is substantial longitudinal curvature in horizontal ice flow, we expect the distribution of stresses to vary locally and influence drainage behavior. The degree of curvature shifts the location at which margin-parallel shear stress decreases toward the inside of the bend, resulting in enhanced tensile stresses and higher velocities on the outside curve of the bend relative to the inside (Echelmeyer and Kamb, 1987). Additional work is necessary to quantify local flow curvature in the vicinity of the water-filled crevasse groups and its relationship to drainage dynamics.

Over the course of the melt season, the configuration of the basal hydrologic network is contingent on the meltwater supply and controls the ice dynamic response to basal lubrication (Weertman, 1958; Röthlisberger, 1972; Walder, 1986; Clarke, 2005; Parizek, 2010; Flowers, 2015). Initially rapid recharge into the subglacial network tends to lead to pressurization and meltwater transport within a distributed or linked-cavity system (Walder, 1986) until turbulent heating develops a more efficient, low-pressure channelized configuration (Röthlisberger, 1972; Clarke, 2005; Schoof, 2010). Furthermore, for short-term variations in water supply comparable to drainage from water-filled crevasses, the ice dynamic response can be significant even when the subglacial network is as efficiently organized as a channelized system (Schoof, 2010; Bartholomew et al., 2012). This is a result of the inability of subglacial conduits to accommodate rapid melt input, requiring an increase in hydraulic gradient, which decreases the effective basal pressure and leads to faster basal sliding (Schoof, 2010; Werder et al., 2013). Therefore, basal effective pressure varies over the melt season as a function of the timing and quantity of recharge. Our model does not explicitly account for the evolution of hydrologic system, but rather assumes a pressurized water layer at the time of surface drainage and basal recharge, consistent with borehole measurements in the region indicating the ice column is near flotation (Lüthi et al., 2002). Observed drainage rates of water-filled crevasses range from 0.001 to 0.012 km2d−1 over approximately 2 weeks to a month (Lampkin et al., 2013). Some water-filled crevasse systems have been observed to drain over shorter time intervals (< day) (Lampkin et al., 2013; Cavanagh et al., 2017; Joseph and Lampkin, 2017). There is considerable inter-annual variability in drainage events over these water-filled crevasse systems. Multiple drainages within a single melt season are correlated with increases in average summer temperatures (Joseph and Lampkin, 2017). During more intense melt seasons, the opportunity for meltwater production increases. This allows for the replacement of infiltrated meltwater within the depressions where ponding occurs. We speculate that both solar insulation and sensible heat flux are the primary factors that drive melt rates. Preferential melting driven by solar insulation as a function of mean azimuth angles within the ponded region may account for differences in melt production across the pond systems, but seasonal variability in melt intensity regulated by sensible heat fluxes may account for fluctuations in inundated area and refilling rates. Drainage from these systems can represent both rapid and protracted meltwater input to the subglacial hydrologic network. Our results also indicate that, in the near future, the areal extent of inundation may be most important in driving regional ice flow response to hydrologic shear weakening, particularly if the observed trend in inundation area documented in this analysis continues (Figure 4B). Although we focused only on meltwater injection via water-filled crevasse drainage, direct surface meltwater infiltration from surface runoff (Figure 1) may be substantial if sustained over the melt season. More work is required to incorporate the evolution of the subglacial hydrologic network and decompose the relative impact of various hydrologic modes of meltwater input to the bed.

Though this analysis has focused mainly on the speed-up of extra-marginal ice due to hydrologic shear weakening, our results demonstrate this process could also lead to a ~7% increase in ice flow within the main trough of Jakobshavn Isbræ (Figure 5A). This impact is small compared to oceanic forcing of the calving front, which drives greater changes in flow speeds near the terminus of the glacier. Though we demonstrate a hydrologically enhanced mass flux of ~16 Gt a−1 through the southern margins (Figure 2A) for a 50% increase in pond extent, it is clearly not enough to account for the total discharge from Jakobshavn estimated at ~ 48 Gt a−1 (Muresan et al., 2016). Bondzio et al. (2017) argue that faster flow due to calving can weaken the shear margins through strain softening, thereby allowing for widespread acceleration outside of the ice stream. In their modeling study, the impact of enhanced fracturing within the shear margins was parameterized using a “damage” factor, which represents the proportion of load-bearing surface area within the margins due to fractures (Borstad et al., 2013; Bondzio et al., 2016). However, stress transmission from an ice stream to adjacent shear margins is known to concentrate into a stress accommodating layer near the base of the ice column where basal drag compensates for enhanced lateral shear stress (Whillans and van der Veen, 2001), suggesting that (Bondzio et al., 2017) damage parameterization would result in an overestimate for the magnitude of extra-marginal ice acceleration that can be attributable to reduced ice viscosity alone. Rather, it is likely that multiple processes are required to fully explain the observations.

Enhance surface melt water infiltration is expected to drive faster summer velocity but would decrease during the winter due to the development of an efficient subglacial drainage system (Sole et al., 2013). This trend may not be the case for regions of the ice sheet above the equilibrium line. Doyle et al. (2014) demonstrate that at higher elevations enhanced summer melt and infiltration are not at levels sufficient to drive the development of efficient subglacial hydrologic system resulting in continued ice flow acceleration into the winter. By 2060, it is expected that supraglacial lakes and ponded features will expand inland by 53% (Leeson et al., 2015). Given the varied conditions in subglacial hydrologic development, the impact of future inland expansion of surface melt water input to the bed via lake drainages will likely result in spatially heterogeneous ice dynamic response across the ice sheet. Under such a scenario, we expect water-filled crevasses to increase in size and distribution throughout Greenland's outlet glaciers, leading to an expanding zone of dynamic influence. This expansion would certainly amplify shear weakening through both basal lubrication and cryo-hydrologic warming from pond drainage and infiltrated surface runoff. We would expect the impact of oceanic forcing of terminus acceleration to provide persistent background strain softening of the near-terminus region, where the direct dynamic impact of calving is most significant (Joughin et al., 2010b), along with upstream accelerations associated with the loss of downstream rigidity. As the terminus retreats, the effects of flow softening may become dominant farther inland. With an increase in surface melt production and runoff, we expect for hydrologic shear weakening due to basal lubrication and runoff infiltration to remain significant farther upstream of the terminus. The spatial and temporal scales over which these various processes operate as well as feedbacks among them are not well understood, requiring further investigation.

Conclusion

In this analysis, we focused only on shear-margin weakening of a glacier flowing through a deeply incised basal trough with well-defined margins. There are other outlet glaciers with a range of boundary and surface mass balance conditions that maintain surface hydrologic systems that deliver meltwater into regions of fast flow. Major marine-terminating outlet glaciers throughout the ice sheet have water-filled crevasses and lakes distributed in various configurations. The type of hydrologic system and the relative location of these features can be important parameters in modifying glacial force balance, resulting in either amplifying or abating mass discharge. With warming regional temperatures, we estimate an increase in the ponded area by up to 30% beyond current extents, accompanied by an increase in stored volume and periodicity of drainage into the shear margins. Under this scenario, acceleration of Jakobshavn Isbræ's extra-marginal ice flow could exceed the estimated 0.6 Gt yr−1 of enhanced mass flux into the ice stream. The long-term impacts are constrained by the rate of subglacial hydrologic evolution with increasing meltwater injection, which will likely be spatially heterogeneous. This process is broadly applicable to many other marine-terminating outlet glaciers. The response to direct meltwater drainage into the shear margins of Greenland's outlet glaciers will largely be dependent on the combination of the unique boundary conditions intrinsic to each of these systems and the magnitude and frequency of surface hydrologic forcing.

Understanding the relative partitioning of observed changes in the basin-wide mass flux requires models that more accurately represent processes that can impact ice flow. For Jakobshavn Isbræ in particular, there is a need to develop a comprehensive model that can simulate basal lubrication from crevasse pond drainage, runoff infiltration and terminus-induced strain softening. An improved understanding of how these various processes interact is required to ultimately understanding their collective impact on the spatiotemporal evolution of the overall flow system.

Author Contributions

DL conceived the basic research concept, developed numerical models, performed data processing, analysis, and preparation of manuscript. BP contributed to development of numerical model components and design, participated in analysis, interpretation, and preparation of manuscript. HS and EL provided ISSM technical support and data for development of numerical models. JC and CJ supported data processing and analysis efforts.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported under National Aeronautics and Space Administration grants NNX14AO68G (DL, JC, and CJ), NNX15AH84G (BP). Additionally, work has been supported under U.S. National Science Foundation grants PLR-1443190 (BP), and ANT-0424589, AGS-1338832 (BP). HS was supported by grants from NASA Cryospheric Science and Modeling, Analysis, and Prediction Programs. We would also like to thank Mathieu Morlighem (ISSM development team member) for access to boundary data.