Elastic deformation plays a non-negligible role in Greenland’s outlet glacier flow

Future projections of global mean sea level change are uncertain, partly because of our limited understanding of the dynamics of Greenland’s outlet glaciers. Here we study Nioghalvfjerdsbræ, an outlet glacier of the Northeast Greenland Ice Stream that holds 1.1 m sea-level equivalent of ice. We use GPS observations and numerical modelling to investigate the role of tides as well as the elastic contribution to glacier flow. We find that ocean tides alter the basal lubrication of the glacier up to 10 km inland of the grounding line, and that their influence is best described by a viscoelastic rather than a viscous model. Further inland, sliding is the dominant mechanism of fast glacier motion, and the ice flow induces persistent elastic strain. We conclude that elastic deformation plays a role in glacier flow, particularly in areas of steep topographic changes and fast ice velocities. Ice flow dynamics in Greenland’s outlet glaciers are influenced by elastic deformation, both in the area of tidal influence up to 14 km inland from the grounding line and further upstream, suggest analyses of GPS observations and numerical simulations.

T he Greenland Ice Sheet has experienced increased mass loss due to global warming over the past decades 1 . Warmer air and ocean temperatures continue to destabilise outlet glaciers and enhance ice mass loss 2 . Projecting their future contribution to sea level changes under different climate scenarios has become one of the most fundamental questions for coastal regions [3][4][5] . While much progress has been made in large-scale ice sheet models, the behaviour and accurate modelling of fastflowing ice streams remain a high priority for simulating future ice mass changes 6,7 . The ice discharge of these fast-flowing outlet glaciers into the ocean is a fundamental process in the ice sheet mass-balance 8,9 . Up to now, little is known on how short-term stress changes control outlet glacier dynamics and hence ice discharge. To model the effect of short-term stress changes, an appropriate material model needs to be employed to represent viscous long-term ice flow and elastic short-term fluctuations of the stress distribution. Current ice flow models discard the elastic part in the underlying material model due to its higher computational demand and complexity 10 . The obvious question is to what extent elastic deformation is crucial for glacier motion and in which areas elastic quantities become important.
In fast-flowing ice streams, sliding is a major component that determines the amount of ice discharge 11,12 and an improved understanding of basal processes is needed to obtain better estimations of the ice outflow 13,14 . Our current understanding of basal conditions and their parameterisation in ice sheet models limits the ability to model the evolution of ice sheets over the next century 15 . The subglacial hydrological system is an important factor controlling the stress at the ice base, as a higher water pressure lubricates the base and enhances sliding speeds [16][17][18][19] . For marine terminating glaciers, a connection between the subglacial hydrological system and the ocean water exists. Ocean tides modulate the water pressure at the grounding line resulting in short-term changes of basal stress. These short-term stress changes are likely transmitted inland and impact the ice dynamics upstream 20 . As direct observations at the base are sparse, numerical models are valuable tools to investigate the variation of stresses.
While it is widely accepted that ocean tides induce nonnegligible elastic strains, other potential sources for short-term changes in stresses exist. Most likely steep slopes or steps in bed topography, where ice moving across bed undulations experiences changing stresses, a prerequisite for developing elastic strains. However, it remains unclear whether and where elastic strains exist within the ice for longer times or decline such that the viscous glacier flow prevails. The impact of elastic contributions related to abrupt stress changes on ice dynamics has not been quantified yet, either numerically or through observations. These elastic strains are related to stresses based on the constitutive equation and are commonly used to indicate material failure (e.g., crevasses, damage, cracks). Interestingly, extensive crevasse fields occur in fast-flowing outlet glaciers in Greenland, which exert an important control on ice discharge towards the ocean. Approaches of linear elastic fracture mechanics are promising to represent rapid failure and crevasse formation in largescale ice sheet models 21 . With increasing tensile stresses above a threshold, brittle failure sets in 22 . Based on these considerations, several logical questions arise: Are the crevasse fields located in regions in which elastic strains and stresses do not decline due to stress changes? Could a viscoelastic model provide novel insights into the role of elastic strains and stresses in Greenland's outlet glaciers?
So far, observations of the elastic effects have primarily been focusing on tides: In Antarctica, it has been shown that tides can alter glacier flow up to 80 km upstream of the grounding line of large ice streams 20,23-27 . The tidally-induced vertical motion of ice shelves and floating ice tongues cause changes in the stress distribution and contributes to the variability of glacier flow near the grounding line 26,28 . A case study of Rutford Ice Stream showed that only a combination of a viscoelastic material model and a non-linear friction law describes observed fortnightly glacier speed modulation appropriately 24 . The viscoelastic flowline simulation of Bindschadler Ice Stream, West Antarctica, indicated basal motion consistent with flow dominated by basal sliding over a relatively weak till (determined by higher exponents in sliding laws) and not by internal deformation 25 . These findings suggest a tidal modulation of the glacier motion facilitated by basal sliding. Despite the broad impact of tides on Antarctic ice streams, only a few studies have considered tidally modulated glacier motion on the Greenland Ice Sheet (GrIS) 29,30 . Observed bending of the ice tongue at Jakobshavn Glacier (at that time 8 km long) reflected elastic and viscous deformation, while fatigue failure caused by tidal flexure results in deep basal and surface cracks along the grounding line 29 . Although tidal stresses only affected the glacier terminus at Helheim Glacier, observations after a calving event showed changes at least 12 km upstream of the grounding line 30 . The tidal range around Greenland varies between 0 and 6 m 26 with a magnitude that is similar to the tides around the Antarctic Ice Sheet. Observations at a marine-terminating glacier of the Qaanaaq ice cap revealed semi-diurnal speed peaks coinciding with low tides 31 . Near the calving front of Bowdoin Glacier, tidal forcing and surface speed were in anti-phase. An underprediction of the amplitude of the semi-diurnal variability is presumably related to either inaccuracy in the surface and bedrock topographies or mechanical weakening due to crevassing 32 .
We employ a viscoelastic model underpinned by an accurate radar-derived ice geometry to investigate the elastic contribution along the central flowline of Nioghalvfjerdsbrae (79°NG), northeast Greenland. The parameters for the simplest homogeneous isotropic viscoelastic material are the viscosity, determining the rate-dependent viscous behaviour, and two elastic constants, namely Young's modulus, affecting the stiffness, and Poisson's ratio, specifying the lateral contraction (see "Methods" section). We rely on observed displacements of GPS stations (see "Methods" section) to assess the accuracy of simulated displacements and validate the material rheology. To model basal sliding we have to know the effective water pressure of the subglacial hydrology and an inferred friction coefficient field (see "Methods" section). The effective pressure is the difference between the ice overburden pressure and the subglacial water pressure and results from the confined-unconfined aquifer system (CUAS) model 33 . Inversion of observed surface velocities provides the friction coefficient field. Based on these results, we show that ocean tides change the ice flow only a few kilometres upstream of the grounding line at 79°NG. Crevasse fields further upstream thus require another mechanism that sustains elastic deformations. A simulation without tides enables us to evaluate the elastic contribution to the glacier flow at 79°NG. We show that sharp bed undulations and glacier motion cause elastic deformations, which do not decline over time. From the gained insights focusing on 79°NG, we discuss how this study sheds light on the elastic deformation of other fast-flowing outlet glaciers of the GrIS (viscous model) as the applied viscoelastic Maxwell model demands that elastic and viscous stresses are the same. In the end, we propose ideas for future modelling work.

Results and discussion
Tidal modulation of glacier flow. For the 79°NG, an outlet glacier of the Northeast Greenland Ice Stream (NEGIS), we first analyse the influence of tides on glacier flow (Fig. 1). In the summer of 2017, we deployed four GPS receivers along a central flowline to measure the tidal influence on glacier flow. The instruments were positioned on the floating part of the ice tongue (GPS-shelf), on the~4 km wide hinge zone (GPS-hinge) that experiences bending with tides (see "Methods" section), and some 14 km (GPS-GL-14) and 45 km (GPS-GL-45) upstream of the grounding line (defined as the upper limit of the hinge zone in 2017, Supplementary Fig. 1). The measured flow velocities of the GPS stations positioned along a flowline ( Fig. 1 and Supplementary Data 1) show an increasing trend towards the grounding line from 290 m a −1 at GPS-GL-45 over 950 m a −1 at GPS-GL-14 to 1460 m a −1 at GPS-hinge and slightly decreases to 1310 m a −1 at GPS-shelf 34 . These observed velocities are mean values measured over a 14 day period during the summer that are potentially about 10% above winter velocities 35 , but fit to the simulated surface velocities of the flowline. We recorded three-dimensional surface movements simultaneously over several tidal cycles.
The processed data (see "Methods" section) reveal a tidal signal only in the floating part at GPS-hinge and GPS-shelf (Fig. 2). At 79°NG, we observe no tidal signal at the grounded GPS stations GPS-GL-14 and GPS-Gl-45, which indicates rapid attenuation within a few kilometres upstream of the hinge zone. Though it is not directly comparable, this is in contrast to observations on large ice streams in Antarctica 20,23 where the tidal signal is transmitted tens of kilometres upstream of the grounding line. Inspired by the processes occurring in Antarctica, we deployed the GPS stations and were surprised to see no tidal signal 14 km upstream the grounding line. Besides the tidal signal, a superimposed signal was simultaneously recorded by the GPS receivers ( Fig. 2b-d). We link this to supraglacial lake drainage 35 .
We set up a numerical model and compare the simulated displacements with observations. The momentum balance of this simulation is solved for a viscoelastic Maxwell rheology 36 (labelled COMice-ve, details see "Methods" section). This viscoelastic model allows us to simulate and distinguish shortterm elastic and long-term viscous responses to external loads like ocean tides [37][38][39] . For the floating ice tongue of 79°NG, we account for observed ocean tides of GPS-shelf in the basal boundary condition. Additionally, we include this tidal signal in the hydrological model CUAS.
The simulated glacier response using a viscoelastic material model reveals an excellent agreement with the observations. We choose elastic constants for Young's modulus of E = 1 GPa and Poisson's ratio of ν = 0.325, which are consistent with other studies in Antarctica 37,39 . Without further tuning, the correspondence of model results and observations is high with differences less than 1% in the horizontal and vertical direction over large parts of the considered time interval ( Fig. 2 and Supplementary Data 2). The free-floating ice tongue moves up and down caused by ocean tides and the simulated detrended horizontal and vertical displacement coincide with the observations (Fig. 2d). At both grounded GPS stations GPS-GL-14 and GPS-GL-45 the simulations show no measurable tidal modulation, which agrees with the observations (Fig. 2a, b). We inferred that the high agreement at GPS-hinge solely arises from the elastic contribution in COMice-ve. Horizontal and vertical displacements computed with purely viscous rheology are not capable of reproducing the observed amplitudes or phase responses (Fig. 2c). Although we varied viscosity and basal sliding parameters (Supplementary Discussion 1), those simulation results cannot reproduce the displacements accurately. That means that the bending and modulation of flow speeds in the hinge zone are a result of the viscoelastic nature of ice. Hence, at least a viscoelastic Maxwell material model (two parameter fluid) is required on the lower part of the glacier to match the observations. The question arises of how a viscoelastic material model and tides influence basal sliding. To what extent are these findings (elasticity influences the glacier flow) only applicable to tidally modulated flow? Other numerical studies also discussed the influence of elastic deformations caused by short-term variations. For instance, large calving events in Greenland cause first a reversal of flow and elastically compress the ice front as large icebergs rotate 40 . The reversal of the ice flow typically lasts for a few minutes and is then followed by a flow acceleration, increased longitudinal strain rates and enhanced responses to tidal forcing 30 . Hence, sudden short-term changes like calving events or lake drainage will induce an additional elastic response, which should be investigated in future work and could be the explanation for flow deceleration and acceleration.
Effect of tides on sliding. We intend to get a complete picture between GPS-GL-14 and the grounding line to estimate how far tides reach inland at 79°NG and influence its basal sliding behaviour. Different basal friction laws exist and need to be handled carefully as basal friction is a crucial boundary condition 11,41 . The applied friction law is Budd-like and non-linear 12,42 and relates basal shear stress τ b to basal velocity v b , the friction coefficient field k 2 and the effective pressure N (see "Methods" section, Eq. 2). As we assume the inferred friction coefficient field to be constant in time for our study (see "Methods" section), changes in the basal shear stress originate from tidal changes in effective pressure and Yellow lines represent contours of the basal friction coefficient field k 2 . White lines denote the simulated water pressure difference Δp w (high tide minus low tide) in the hydrological system. The flowline profile (red line) is used for a flowline model of glacier motion. The grounding line, where the ice becomes afloat is at 0 km. The hinge zone, the zone in which the ice tongue deforms due to bending, is shown by two black lines, one at its upper limit directly at the grounding line and one at the lower limit 4 km downstream the grounding line. The inset shows the CUAS domain that encompasses the NEGIS catchment including the 79°NG catchment as a solid black polygon. The red rectangle marks the area close to the grounding line shown in the figure. basal velocities. To disentangle the role of each component in the friction law, we analyse the response of the simulated subglacial hydrology to tides resulting in time-dependent effective pressure fields (see "Methods" section). The simulated vertical movement of the floating base caused by tides result in time-dependent displacement and hence velocity changes.
With this tidally modulated setting along the flowline, we show a comprehensive analysis of the glacier system for all quantities included in the basal boundary condition ( Fig. 3 and Supplementary Data 3 and 4). Please note that results are shown point-wise relative to a control run (CTRL) without tidal forcing in the subglacial hydrology system (grounded part) and the water pressure (floating part). The simulated head along the flowline (the corresponding elevation of the water level in a hypothetical borehole connecting glacier surface to bed) is very smooth (Fig. 3b) and spatial variations of the effective pressure primarily originate from the ice overburden pressure and thus ice thickness variations. In the hinge zone (from 0 to 5 km), the simulation reveals a tidal influence on the basal velocity where the highest differences to the CTRL run occur at the transition between low to high tide. The amplitude in the basal velocity behaves similar to the displacements ( Fig. 2 and Supplementary Fig. 2), where smaller amplitudes appear for the purely viscous or a stiffer material (dashed and dotted lines in Fig. 3d). Simulated velocities and displacements show a tidal dependency that drops from ±10% at the grounding line to ±0.2% at GPS-GL-14 ( Fig. 3d) to none for GPS-GL-45 (not shown). For the grounded part, the transition from low to high tides and vice versa induces the largest changes compared to the spatial evolution for all other tidal states.
Why occur the lowest or highest basal velocities not at high or low tide but the transition points ( Fig. 3d)? Although we expect the tidal effects on the hydrological system to be relevant for glacier flow, we conduct a further simulation without tidal modulation to disentangle the effect of tides on the friction law and the rheology. Hence, the effective pressure in the subglacial water system is kept constant (CTRL run) and only the basal boundary condition at the floating tongue experiences tidal modulation. The velocity variations induced by the vertical tidal movement of the floating tongue reaches around 10 km upstream the grounding line. The bending of the floating tongue influences the basal velocity and high tide leads to lower velocities at the base of the grounded ice, while low tide causes higher velocities ( Supplementary Fig. 3). This scenario highlights that the tidally dependent change of the hydrological system alters the flow behaviour: the envelopes of the velocity changes occur for the transition between low and high tides (Fig. 3d). If we include tides in the hydrological system, the effective pressure reacts with a decrease in magnitude and a phase shift to tidal forcing depending on the distance to the grounding line (Supplementary Discussion 2). We find the effect of tidal forcing on the hydrological system is limited to a small band up to 20 km upstream the grounding line where the amplitude of the tidal signal decays to only 10% of the signal at the grounding line (Supplementary Figs. 4 and 5). Caused by tidal forcing, relative changes in the effective pressure are small, except at the grounding line ( Fig. 3e and Supplementary Figs. 4 and 6). Additionally, tidal-related changes in the simulated basal velocity are only visible over a distance of less than 14 km upstream the grounding line (Fig. 3d), which is consistent with the observations at GPS-GL-14 (Fig. 2). To validate a non-linear friction law at 79°NG similar to the previous work done in Antarctica 37 , we have to consider grounded GPS stations closer to the grounding line that record tidal modulations.
Previous observations of radio-echo sounding data from NEGIS suggested that upstream areas of the ice stream are controlled by variations in basal lubrication while downstream regions are confined by basal topography 43 . The observation of the downstream region fits the simulation results of the central flowline model for 79°NG shown here. In the basal shear stress, we see impacts of tides with changes up to 40% to the CTRL run near the steep increase in the bed slope at −9 km, where also the subglacial hydrology reveals tidal effects (Fig. 3c, e). In the last ~5 km upstream the grounding line, the friction coefficient is small and damps changes in the basal shear stress due to tides (Fig. 3b). Low friction coefficients facilitate the fast response of glacier speeds to changes in the basal water pressure, which is in line with areas of seasonal speedup and short pulses of acceleration arising from supraglacial meltwater input modulating basal frictional stresses 35 . The tidal variation of the hydrology causes a superposed variation in the basal velocity and consequently a larger change in basal shear stress than neglecting ocean tides in the hydrological system. Overall, ocean tides alter the glacier flow in our case not as far inland as in Antarctica 20 . In our simulations, the area of high transmissivity at 79°NG ( Supplementary Fig. 6b) is confined by the steep basal topography −9 km upstream of the grounding line. Large gradients in the bed or surface topography lead to large gradients in ice overburden and effective pressure. These increase channel wall melt rates 17 and therefore transmissivity 33 . The simulated low aquifer transmissivity at~10 km and further upstream the grounding line hinders the pressure wave to travel further inland at 79°NG (Supplementary Discussion 2). This matches very well with our measurements at GPS-GL-14 where no tidal signal was observed.
What governs glacier motion. Ice sheet models typically neglect the elastic deformation and compute ice sheet flow based on a viscous model. To estimate elastic contributions to glacier flow, it is necessary to get better insights into processes influencing glacier motion. We analyse the role of sliding and vertical shear deformation independent of tides for 79°NG. Based on this analysis, we aim to derive implications for the entire ice sheet of Greenland. The analysis in this section uses (i) a viscoelastic simulation along the flowline of 79°NG and (ii) a continental scale ice sheet simulation. For the latter, we rely on a viscous model since a viscoelastic simulation is unavailable. The overarching goal is to identify regions where the consideration of elasticity is recommended for a better understanding of outlet glacier dynamics. We analyse vertical profiles of velocity ratios determined by the velocity v = |v| at each point relative to its basal velocity v b = |v b | (Fig. 4a). The simulated velocity field differs from a sliding dominated plug flow (grey area) in regions where the contribution of vertical shear deformation is enhanced. Topographic bed changes seem to induce these higher contributions (e.g., at distances of −18 and −9 km). Slightly before the bed bump peak at −9 km, simulated vertical shear deformation is highest and topographic changes slow down the ice at the base (velocity at the surface is larger than at the base). At this location, we find a maximum of 25% higher velocities at the surface. At some locations, the velocity is even lower at the surface than at the base. Such a feature is known from viscoelastic flow over an undulated bed 44,45 . The surface velocity is up to 5% lower than the basal velocity. In general, plug flow is the dominant mechanism and vertical shear deformation is small. The simulation results fit radar and seismic observations indicating deformable subglacial sediment in the far upstream part of NEGIS, which facilitates sliding 46,47 . At 79°NG, the simulation results support that fastmoving outlet glaciers flow mainly by basal sliding 30 .
To get better insights into the role of deformation in glacier motion, we first analyse the ratio of vertical shear deformation (d) to sliding (s) and second the amount of viscous to total strain ( Fig. 4b and Supplementary Data 5). We introduce a measure M a d=s 2 ½0; 1 representing the absolute contribution of vertical shear deformation to sliding with the ice thickness H, the surface and bed elevations h s , h b , respectively. This depth-integrated measure is zero for a sliding dominated plug flow and allows us to identify regions in which shear deformation of the vertical ice column is pronounced. Overall, the magnitude of M a d=s is higher near steeper bed slopes (Fig. 4b). Vertical shear deformation does not exceed 22% along the entire flow line, although this profile contains a rather steep basal topography. The visualisation of this measure shows a horizontal shift to the downstream side of the slopes (Fig. 4b).
This offset occurs as vertical shear deformation is highest at the end of steep topographic changes indicated by declining slope values. If the friction coefficient reaches low values (near the grounding line or at a distance of −35 km), vertical shear deformation to sliding M a d=s is small (Fig. 4b). The reduction of friction facilitates plug-flow. In the remaining regions, a deformation dependency on slope and friction coefficient is apparent. For a viscous simulation of GrIS motion 48 , we also consider the contribution of vertical shear deformation to sliding (Fig. 5a). Sliding is dominant in all fast-moving areas but also in the interior of Greenland. Vertical shear deformation occurs up to 35% in the interior or at thin margins where ice is presumably frozen to ground 49 .
The remaining question is whether bed undulations can induce elastic responses that persist and consequently contribute to glacier motion. The viscoelastic Maxwell model for the flowline of 79°NG allows us to distinguish between the viscous and elastic strain contribution as both sum up to the total strain 36 . We exploit this property and quantify the ratio of viscous to total strain. We define the number M v e=v 2 ½0; 1 based on Euclidean matrix norms of the viscous (ϵ v ) and elastic (ϵ e ) strain This measure M v e=v reaches the value 1 for vanishing elastic strains. Without short-term changes of tides, the amount of viscous relative to total strain consists of 6% elastic strain on average to a maximum above 30% for the grounded part (Fig. 4b difference of M v e=v to 1, Supplementary Fig. 7). The ratio of viscous to total strain is close to 1 at the free-floating ice tongue, meaning that the viscous flow dominates as no short-term changes happen. Once elastic strains caused by transient initialisation vanish, the partition of viscous to total strain M v e=v remains similar for long simulation times (Fig. 4b for 1a, Supplementary Fig. 7 for up to 10a). However, the ice flow over undulated beds induces unexplored time-independent elastic deformations occurring at grounded parts of the fastflowing 79°NG.
So far, large ice sheet models did not consider elastic contributions. The implication to neglect this component is not clear yet. For a glacier motion above undulated beds, the elastic strain increases and kinematics lead to higher displacements (integral of the strain), which enhance the flow velocity of the glacier (Supplementary Fig. 8). The results of the viscoelastic flowline simulation reveal the existence of elastic strain for changing bed slopes (Fig. 4b). At 79°NG, bed slopes are too fluctuating for the elastic strain to vanish (which is due to Blueish colours depict velocities that are more than 2% smaller than its basal velocities, while reddish colours show domains in which v is at least 2% higher than v b . The blue lines (right axis) depict the simulated surface velocities compared to MEaSUREs velocity 64,65 . b The reddish lines show the normalised friction coefficient and effective pressure along the grounded part of the flowline while the blueish curves present the contribution of viscous to total strain M v e=v and the measure of vertical deformation to sliding M a d=s that reacts with a slight horizontal delay to the normalised slope of the base shown in grey. All those quantities correspond to one year of simulation time.
fluctuating stresses), except near −4 or −40 km where M v e=v is almost 1. Consequently, a rougher bed fosters more elastic deformation. For a viscoelastic Maxwell model, the total strain is the sum of elastic and viscous strain (see "Methods" section). The elastic strain is around one order of magnitude smaller than the viscous strain (Supplementary Methods 1 and Supplementary  Fig. 9). For instance, a constant elastic strain of 2‰ over 50 km induces an additional elastic elongation of 100 m. In the end, a viscoelastic simulation reveals a higher velocity for a fast-flowing outlet glacier than a viscous one caused by the additional elastic deformation.
We aim to transfer local insights at 79°NG to Greenland by determining regions in which elastic deformations may be crucial for understanding outlet glacier dynamics. To deal with the lack of a viscoelastic simulation for entire ice sheets, we leverage the fundamental character of a Maxwell rheology: elastic and viscous stresses are equal. Simulated stresses obtained by a viscous model (initialisation experiment 48 ) can be used to compute principal (i.e., largest) elastic strains (see "Methods" section) on the entire surface of the GrIS (Fig. 5b-e and Supplementary Fig. 10). Without short-term fluctuations, we would expect that the elastic strain declines with time and vanishes. The largest elastic strain is found in fast-moving outlet glaciers and reaches at most 1.14‰. As the elastic strain does not vanish for Greenland (Fig. 5b-e and Supplementary Fig. 10), we analyse the change in first principal stress with respect to time (see "Methods" section). Changing first principal stress over time often coincide with regions of large elastic strains in the fast-moving parts of outlet glaciers ( Fig. 6 and Supplementary Fig. 10). In these areas, elastic strains did not decline as changing stresses prevent this. The estimated occurrence of elastic strains at 79°NG (Fig. 5c) matches with the spatial extent of time-independent elastic strain from the viscoelastic flowline simulation (Fig. 4b). This agreement indicates that slope change and viscous flow above 250 m a −1 cause stress changes inducing elastic strains. For assessing the area of Greenland's ice affected by an elastic contribution arising from bed undulations, we use the BedMachine bedrock topography 50 and compute the directional slope ( Supplementary  Fig. 11). A slope larger than 5% is present in around 10% of the grounded ice of GrIS. For the flowline simulation of 79°NG, the mean slope of the last 45 km upstream the grounding line is 4.7% where at least in this region the elastic strain is not negligible (Fig. 4b).
How does this time-independent elastic strain compare in magnitude to tidal forcing? We first estimate the stress change caused by ocean tides ranging from 0 to 6 m around Greenland (Fig. 6). For the three floating tongues glaciers, tides span from~0.5 m for Ryder Glacier,~1.0 m for 79°NG tõ 2.2 m for Petermann Glacier. Due to hydrostatic equilibrium, a diurnal tidal range induce 5-20 kPa d −1 stress change at the floating tongues. At 79°NG, the estimated stress change is 10 kPa d −1 = 3650 kPa a −1 at the glacier tongue. Independent of tides, we reach one third of this stress change magnitude at the grounded part of outlet glaciers for Jakobshavn Isbrae, 15% at Kangerlussuaq Glacier and 10% at Helheim Glacier (Fig. 6). We find in other outlet glaciers close to the grounding line about 10% of this value declining to 1% over some tens of kilometres. In the flowline simulation of 79°NG, elastic strains occur in sliding dominated flow and in regions with higher vertical deformation (Fig. 4). In the end, faster flow over rougher beds induces higher elastic strains caused by higher stress changes. Caused by the motion of fast-flowing outlet glaciers over undulated beds, the stress changes are lower of those due to tides. While tides force an elastic elongation and compression on daily time scales, the motion of fast-flowing outlet glaciers induces an elastic elongation independent of time.
As there is no way to conduct in-situ observations of stress, we relate elastic deformation to the occurrence of cracks. Cracks are a representation of the solid nature of a material. Brittle failure is determined by a small process zone in which local elastic quantities causes a fracture, not viscous ones 51 . The extensive crevasse fields that we observe in outlet glaciers in Greenland are illustrative proofs of this elastic contribution ( Supplementary  Figs. 10 and 12). Furthermore, crevasse fields in ablation zones facilitate surface meltwater discharge to the ice base. This contributes to seasonal speed-up via lubrication, which increases elastic deformation. Damage in the lower tens kilometres of tidewater glaciers also controls calving dynamics of tidewater glaciers and thus ice discharge rates as feedback 22 . Taking into account that Jakobshavn Isbrae and Kangerlussuaq Glacier contributed to an estimated 8 mm of global sea-level rise over the last century 2 , understanding their future response is crucial. All ice discharge of the GrIS and in particular Jakobshavn Isbrae and Kangerlussuaq Glacier (Fig. 5b, d and Supplementary Fig. 13) takes place in a settings where the bed is lubricated 52 , sliding speeds are high (Fig. 5a) [53][54][55][56][57] , surfaces are crevassed (Fig. 5b- Supplementary Fig. 12).
Crevasses also occur in slow-moving outlet glaciers in which the basal topography will lead to necessary stress changes (Fig. 6) and to time-independent elastic strains ( Supplementary Fig. 10). The extensive match of massive crevasse fields with areas of large stress change supports the claim that glacier flow over undulated beds induces persistent (time-independent) elastic strains.
Conclusion. By comparing GPS observations with numerical modelling, we demonstrated the need to describe tidal effects on ice by a viscoelastic Maxwell model to reproduce glacier dynamics in the lower part of the 79°NG. The tidal influence, not captured by a viscous model, is limited to the hinge (bending) zone downstream the grounding line and up to 10 km upstream the grounding line. For the grounded ice, the subglacial hydrological system is relevant as it causes a phase delay in basal velocities. In the absence of tides, we found sliding to be the dominant mechanism of fast-moving glacier motion, and a absolute contribution of vertical shear deformation to sliding of less than 22% for the viscoelastic simulation of 79°NG up to 35% for the viscous simulation of Greenland. Deducing from our viscoelastic modelling, the ice flow induces a persistent elastic to total strain of 6% on average up to a maximum of 30% for the considered flowline of 79°NG.
Based on the viscoelastic simulation of 79°NG, we developed a new perspective on deformation in outlet glaciers by transferring our findings from the flowline model to a continental scale simulation. We demonstrated that elastic deformation additionally contributes to the motion of outlet glaciers. Independent of tides, elastic deformation plays a role in areas of non-steady stress fields, such as near steep changes in topography or as ice flow speeds up. These new insights are a step forward to better understand the motion of glaciers since the identified regions are currently poorly constrained in ice flow models. In future works, improving the fidelity of ice flow models by considering the elastic contribution may capture the behaviour of fast-moving outlet glaciers. Under the overarching goal of narrowing uncertainties in sea-level projections, it is arguably important to assess how ice discharge to the ocean is affected by including an elastic deformation in the overall ice flow.

Methods
Viscoelastic modelling. The simulations COMice-ve were conducted with the finite element programme COMSOL applying a viscoelastic Maxwell material model that solves for the unknown displacements (along-flow and in thickness direction) 36 . We prescribe the ice geometry along the flowline from a highresolution basal topography obtained from airborne ultra wide-band radar measurements. All simulation results are obtained using this high-resolution geometry subject to a non-linear Budd-like friction law 42 . The position of the grounding line is the lowest position of the landward limit of the hinge zone defined from SAR interferometry. We show that this landward limit (upper flexure limit) does not move more than 500 m in the time period of observations. Such small changes make no difference in the simulation results, and for simplicity we assume a fixed grounding line at the lowest position of the upper flexure limit derived by SAR interferometry.
COMice-ve solves the underlying viscoelastic Maxwell equations in a small strain setting that is valid for simulation times up to few years. We discretise the model domain with a triangular mesh with a horizontal resolution of 100 m. This fine spatial resolution is necessary to get rid of the mesh discretisation error ( Supplementary Fig. 14). The shape functions for the unknown displacements are quadratic Lagrange polynomials. The viscous strain components are additional internal variables for the Maxwell rheology and we use shape functions of linear discontinuous Lagrange type to save computational effort. The viscous and elastic stress of a Maxwell material are the same while the strain is the sum of elastic and viscous strain. The material relation of the volume-preserving deviatoric stress (σ D ¼ σ À 1=3 trðσÞI) to elastic strain ε D e or viscous strain rate _ ε D v is given by the constitutive equation with the elastic material constant μ = E/(1 + ν) 10 . The viscosity η is non-linear and based on Glen's flow law using a temperature field from previous modelling 61 . We choose the same model parameters as used in previous work 36 . The major difference to previous applications investigating calving of ice shelves 36,62 is the incorporation of a friction law. The general friction power law (Budd-like) relates the basal shear stress (τ b ) to the sliding velocity v b with the stress exponent m. This exponent is positive and unequal to 1 for a nonlinear realisation of the friction law. We choose m = 3 and adjust the unknown friction coefficient field k 2 (see section "Inversion for basal friction coefficient" in "Methods" section) to result in a modelled flow velocity that best matches the observed surface velocity 37,63 . For the effective pressure N = p i − p w , the difference of ice overburden pressure p i and the subglacial water pressure p w , and the friction coefficient field k 2 we rely on external products. The effective pressure is obtained from a subglacial hydrology model. The output is used to perform a one-way coupling by feeding the resulting spatial and temporal distribution of the effective pressure of the subglacial system along the flowline back into the friction law (basal boundary condition for grounded ice in the viscoelastic model). The friction coefficient field is inferred by solving an inverse problem for the nonlinear friction law. The inferred friction coefficient k 2 inv is related to a simpler effective pressure parameterisation N inv applied for the inversion. Consequently, we adjust the unknown k 2 to include in the friction law of the viscoelastic model by the following rule: This re-scaling ensures that the basal shear stress fed to COMice-ve is comparable to the basal shear stress, which is derived from inversion. The boundary condition at the ice-atmosphere interface is traction-free and at the inflow boundary, we take a constant velocity based on the MEaSUREs velocity data-set 64,65 . Given the observed vertical displacements of the GPS-shelf, we can compute the tidal pressure field perturbation at the base of the ice tongue. The floating ice mass is subject to a hydrostatic water pressure field normal to the submerged ice-ocean boundary line, depending on tide height (e.g., high tide leads to larger water pressure). We put much effort in reproducing the observed velocities in the grounded part simulated by COMice-ve. However, the 2D flowline model lacks lateral resistance exerted by the fjord side-walls at the floating tongue. As this effect is inherently not covered we rely on a buttressing formulation 66 . The lateral resistance is accounted for over the whole floating tongue area (i.e., x > x gl ) by adding a body force f in the momentum balance equation, such that with v x the velocity in horizontal direction. In the grounded part, the lateral resistance coefficient K is zero, while on the floating tongue K > 0. The exponent m lr is set equal to m. Here, we tune the lateral resistance coefficient spatially to obtain a reasonable decline of flow velocities downstream from the grounding line. We conducted the benchmark experiments ISMIP-HOM B and D for verification of the long-term viscous response of the applied viscoelastic model (Supplementary Methods 3 and Supplementary Figs. 8 and 15). The simulated short-term elastic response with the applied model was already verified 62 .
Global positioning system (GPS) processing. The GPS data were processed using the GIPSY-OASIS software package including high-precision kinematic data processing methods 67 with ambiguity resolution using Jet Propulsion Laboratory (JPL)'s orbit and clock products, constraint on kinematic position solution. We use the GIPSY-OASIS version 6.4 developed at JPL, and released in January 2020 68 . We use JPL final orbit products which include satellite orbits, satellite clock parameters and Earth orientation parameters. The orbit products take the satellite antenna phase centre offsets into account. The atmospheric delay parameters are modelled using the Vienna Mapping Function 1 (VMF1) with VMF1grid nominals 69 . Corrections are applied to remove the solid Earth tide and ocean tidal loading. The amplitudes and phases of the main ocean tidal loading terms are calculated using the Automatic Loading Provider (http://holt.oso.chalmers.se/loading/) applied to the FES2014b 70 ocean tide model including correction for centre of mass motion of the Earth due to the ocean tides. The site coordinates are computed in the IGS14 frame 71 . We convert the Cartesian coordinates at 5 min intervals to local up, north and east for each GPS site monitored at the surface of the 79°NG. In addition, we use Waypoint GravNav 8.8 processing software. We applied kinematic precise point positioning (PPP) processing using precise satellite orbits and clocks. The site coordinates are computed in the IGS14 frame and converted to WGS84 during data export at 15 s interval. To avoid jumps between daily solutions of the Waypoint PPP product, as the data is recorded in daily files, we merged three successive files prior to processing to enable full day overlaps. In a second step, the 3-day solutions are combined using relative point to point distances. To avoid edge effects, we combined the files in the middle of each 1-day overlap and removed outliers. To estimate the vertical and horizontal displacements we detrended the vertical, northing and easting components of the data after mapping into a polarstereographic projection and alignment in the ice flow direction. Finally, the detrended data was low pass filtered to suppress signal frequencies larger than 1/3600 Hz and re-sampled to 5 min interval to match the GIPSY-OASIS product. A comparison of both processing software products revealed very similar results, with slightly less noise observed in the GIPSY-OASIS product (Supplementary Figs. 16 and 17). The PPP processing with the GIPSY-OASIS software Package of the GPSshelf data was used as input for the simulation.
Modelling of subglacial hydrology. We applied CUAS to compute the effective pressure N = p i − p w at the ice base. The model uses an equivalent porous medium (EPM) approach and accounts for efficient and inefficient drainage. A Darcy-type groundwater flow equation for the EPM is used to evolve the hydraulic head h, where the transmissivity T locally adjusts based on channel and cavity evolution. Areas of high transmissivity represent efficient (channelised) flow, whereas low transmissivity is interpreted as inefficient (distributed) flow.
We used a new version of the BedMachine geometry for Greenland 50 including the latest AWI airborne data 72 . In addition to ice sheet basal melt from the Parallel Ice Sheet Model (PISM) output 73 , we updated and extended the domain outline. We have interpolated the data sets onto a 1.2 km resolution grid (coarse) for the entire NEGIS area and a 300 m resolution grid (fine) covering the area around the 79°NG (Supplementary Fig. 18). We manually set the domain outlines guided by the BedMachine ice mask, the pressure-adjusted basal temperature contour of 0.1°C of the PISM output, digitised ice front and grounding line positions, and ice flow velocities (≥10 m a −1 ). We used no-flux boundary conditions at lateral margins (homogeneous Neumann boundary conditions for the head) and Dirichlet boundary conditions (zero head) at grounding lines. At the ocean boundary the transmissivity was set to T max ¼ 100 m 2 s À1 while the initial value was 0.2 m 2 s −1 elsewhere. We chose the subglacial water pressure to be equal to the ice overburden pressure as initial condition for the head. We choose the same model parameters as used in previous work 33 .
The output of a 50 years long steady-state run on the coarse grid ( Supplementary Fig. 18) was used to run a 1 year long simulation on the fine grid in the 79°NG nested sub-domain to let the model adjust to the changes in grid resolution. The effective pressure at the end of this spin-up was used to adjust the friction coefficient (Eq. 3). We run the model for additional 19 days with and without tidal forcing at the grounding line with hourly output to provide time series of N along the profile shown in Fig. 1 as input for COMice-ve (Eq. 2). We applied the tidal forcing as a time dependent Dirichlet boundary condition for the head at the 79°NG grounding line.
The simulated subglacial drainage system ( Supplementary Fig. 18a, b) shows relatively low effective pressure and transmissivity at the land terminating ice margin, indicating less efficient drainage. Observations indicate that efficient channels can form in the ablation area at the ice sheet margin 74 . At the margin, the ice is usually thin (decreased creep-closure rates), the surface slope is steep and thus, the gradient of the hydraulic potential is large (enhanced channel wall melt). A sufficient amount of water from surface runoff or lake drainage is available at the glacier base via surfaceto-bed hydrological connections (e.g., crevasses or moulins) to allow efficient channel formation. In the model, only ice sheet basal melt is considered as the supply for the hydrological system 33 . At the land terminating margin, basal melt from the PISM simulation is low or absent and thus efficient drainage is inhibited in our simulations. The transmissivity is high in the vicinity of the grounding line indicating efficient (channelised) water transport and declines with distance upstream. In this area, N is close to zero and increases towards GPS-GL-14.
Inversion for basal friction coefficient. The viscoelastic NEGIS ice flow model setup makes use of a basal friction coefficient field k 2 inv that is retrieved by an inversion method. For the inversion of the basal friction coefficient, we operate the Ice-Sheet and Sea-level System Model 75,76 , an open source finite element flow model appropriate for continental-scale and outlet glacier applications. The modelling domain is the flowline profile (Fig. 4) and the ice flow is modelled by the full-Stokes equations. Model calculations are performed on a structured finite element grid with a horizontal resolution along-flow of 0.2 km. The domain is vertically extruded with 15 layers refined to the base. A thermo-mechanical coupling is not performed, but a realistic ice rigidity is prescribed based on thermal spin-up 77 . We performed the inversion for the two different bed geometries: (1) BedMachine geometry for Greenland 50 including the latest AWI airborne data 72 and (2) the high-resolution basal topography from airborne ultra wide-band radar measurements.
Within the inverse problem a cost function (J), that measures the misfit between observed, v obs , and modelled horizontal velocities, v x , is minimised. We use the observed velocities from the MEaSUREs project 64,65 as target within the inversion. This data-set composes a long-term mean (20a) from several remote sensing products and is independent of our GPS measurements. Please note that v obs is the velocity in the curvilinear coordinate system along the flowline. The cost function is composed of two terms which fit the velocities in fast-and slow-moving areas. A third term is a Tikhonov regularisation to avoid oscillations due to over fitting. The cost function is defined as follows: where ε is a minimum velocity used to avoid singularities and Γ s and Γ b are the ice surface and ice base, respectively. An L-curve analysis was performed to pick the Tikhonov parameter γ t (Supplementary Fig. 19a). We obtained a good agreement to the observed velocities by choosing γ 1 = 10, γ 2 = 1 and γ t = 1 × 10 −8 . The inverse problem is run for the non-linear friction law (Eq. (2)), but uses a simple parameterisation of the effective pressure 78 with N ¼ ϱ i g h þ minð0; ϱ w g z b Þ, where h is the ice thickness, h b the glacier base and ϱ i = 910 kg m −3 , ϱ w = 1000 kg m −3 the densities for ice and fresh water, respectively. The parameterisation accounts for full water-pressure support from the ocean wherever the ice sheet base is below sea-level, even far into its interior where such a drainage system may not exist. With this parameterisation of the effective pressure, the inferred friction coefficient field is independent in time. The obtained results for the inferred friction parameters k 2 inv and simulated surface velocities are shown in Supplementary Fig. 19b, c.
Hinge zone by means of interferometry. We applied Synthetic Aperture Radar Interferometry (InSAR) on Sentinel-1 TOPS SAR data to detect the upper and lower limit of tidal flexure. Interferograms are formed from interferometric wide swath single look complex data acquired at times t 1 and t 2 separated by a temporal baseline of 6 days. Assuming constant horizontal ice flow within a 6 day time period we subtracted another interferogram with data acquired at times t 2 and t 3 to isolate vertical displacements due to ocean tides 61,79,80 . This way we were able to produce 26 double differential interferograms from which we manually delineated the upper and lower flexure limit (Supplementary Fig. 1). To minimise the effect of inter annual grounding line migration when compared to the model results we focus on the year 2017. Along the profile used for the viscoelastic modelling we found a maximum variation of 520 m and 880 m for the upper and lower flexure limit, respectively. These variations result from two main reasons (1) the timing between satellite pass and tide level is always different and (2) the horizontal ice flow is not constant within 6 days and hence residua remain in the double differential interferograms 81 .
Changes of first principal stress for Greenland. The elastic deformation of a Maxwell material is non-negligible when subject to short-term changes in stress, such as those caused by calving events or ocean tides. As a Lagrangian ice parcel passes over an undulating bed, the stress experienced might, however, also change enough to induce non-negligible elastic deformation. Caused by such stress changes (Fig. 6), we considered areas where time-independent elastic strains might occur caused by the viscous glacier flow (Supplementary Fig. 10). Elastic stresses are proportional to elastic strains and viscous stresses are related to viscous strain rates based on their constitutive equations. The elastic strain exists wherever viscous strain rate is present-hence in any regions of non-zero velocity gradients. To sustain elastic strain, the temporal change of stresses is necessary as the long-term stress response of a Maxwell model corresponds to the one of a viscous fluid (Supplementary Fig. 8).
Based on the MEaSUREs velocity field 64,65 we generate a set of more than 60,000 flowlines with a point to point spacing of 200 m. The seed points of the flowlines were distributed every 5 km across Greenland, where the surface elevation exceeded 500 m. Along each flowline, we estimate the change of first principal stress over time using a raster with 500 m pixel resolution interpolated from the simulated first principal stress data set. We choose the pixel spacing to preserve the finer mesh used in the simulation in areas of fast flow. In total, we interpolate the stress change of more than 100 million points back to the initial input stress field using the same inverse distance weighting interpolation scheme.
To estimate the change in water pressure over a tidal cycle around Greenland, we used the Arctic Ocean Tidal Inverse Model 82,83 . We apply the model to a 5 km grid covering Greenland's coastline (including a buffer zone of 50 km towards the ocean) in a temporal resolution of 10 min for a tidal cycle (one month). Finally, the tidal range is given by the difference of the tide level extrema within each Pixel. This allows us to get a rough idea of the change in water pressure and its regional variation around Greenland (Supplementary Fig. 10).