Complex Greenland outlet glacier flow captured

The Greenland Ice Sheet is losing mass at an accelerating rate due to increased surface melt and flow acceleration in outlet glaciers. Quantifying future dynamic contributions to sea level requires accurate portrayal of outlet glaciers in ice sheet simulations, but to date poor knowledge of subglacial topography and limited model resolution have prevented reproduction of complex spatial patterns of outlet flow. Here we combine a high-resolution ice-sheet model coupled to uniformly applied models of subglacial hydrology and basal sliding, and a new subglacial topography data set to simulate the flow of the Greenland Ice Sheet. Flow patterns of many outlet glaciers are well captured, illustrating fundamental commonalities in outlet glacier flow and highlighting the importance of efforts to map subglacial topography. Success in reproducing present day flow patterns shows the potential for prognostic modelling of ice sheets without the need for spatially varying parameters with uncertain time evolution.

H igh spatial variability in the flow of the Greenland Ice Sheet is apparent from observations 1 , and capturing this variability is essential to any modelling effort targeting the future evolution of the Greenland Ice Sheet 2 , yet projections of ice discharge into the ocean remains a major wild card for twenty first century sea-level projections 3,4 . Ice flow is nonlinearly related to ice thickness, with small uncertainties in ice thickness potentially leading to large biases in discharge estimates 5 . Consequently, our ability to reproduce observed flow patterns, and thus ice discharge, is contingent on accurate knowledge of ice thickness; however, limited knowledge of thickness has hindered modelling efforts to date. In Greenland, the NASA airborne mission Operation IceBridge (OIB) 6 has added many 1,000 km of radar-derived ice thickness profiles since 2009, nearly doubling the coverage available at that time. To help fill in remaining gaps, mass-conserving interpolation methods 7 have been used to derive flow-compatible, high-resolution maps of ice thickness and subglacial topography 8 . The result has been a substantial improvement in our knowledge of subglacial topography, particularly in the deep channel-feeding outlet glaciers. At the same time, code parallelization, combined with high-performance computing, have begun to make high-resolution ice sheet modelling tractable.
Combining these advances allows us to pursue a set of numerical experiments to investigate whether spatially complex flow patterns in outlet glaciers can now be captured in whole-icesheet simulations using only ice-sheet-wide (spatially uniform) parameters, without local 'tuning' applied to individual grid cells.
Here we use the Parallel Ice Sheet Model (PISM) 9 , a computationally efficient, three-dimensional model 10,11 , coupled with models of subglacial hydrology and basal sliding, to simulate the velocity field of the Greenland Ice Sheet at high resolution (o1 km). We demonstrate that outlet glacier flow can be captured with high fidelity if ice thickness is well constrained and vertical shearing as well as membrane stresses are included in the model (without solving the full-stress configuration), while computing flow from vertical shearing alone and/or using low-resolution ice thickness leads to poor agreement with observations. Overall root mean squared (r.m.s.) velocity differences decrease with increasing model resolution. This indicates that ongoing improvements in the mapping of subglacial topography, together with improvements in modelling resolution, go a long ways towards improved wholeice-sheet numerical simulations.

Results
Calibration and validation. Using the ice sheet geometry given by ice thickness and subglacial topography from Morlighem et al. 8 (data set MO2014), we calculated diagnostic velocity fields with PISM. We calibrated the model at a horizontal grid resolution of 1,500 m by using 15 runs to explore the parameter space defined by three spatially uniform parameters controlling ice dynamics and basal processes, selecting the run with the best fit to observed outlet glacier flow (see Methods and  Supplementary Table 1).
We assessed model performance by comparing observed and simulated surface velocities along cross-flow profiles of 29 large (flux ZB3 Gt yr À 1 ) outlet glaciers (Fig. 1a, Supplementary Table 2), which account for B2/3 of the ice-sheet discharge (B520 Gt per year in 2008) 4 . Along these profiles, we sampled remotely sensed 12 and modelled surface velocities every 250 m using bilinear interpolation. We quantified the model's ability to capture the spatial variation in flow structure using the Pearson r correlation coefficient and the r.m.s. difference between observations and simulations along the 29 profiles. The correlation coefficient provides an indication of the model's ability to capture the velocity variation along a profile (related to the channelization of flow), while the r.m.s. value captures the ability of the model to match the magnitude of the velocity, indicating how well the model reproduces discharge flux of each glacier.
Within the tested parameter space, 9 of 14 calibration experiments are within 25% of the smallest total r.m.s. difference (that is, the r.m.s. difference between observed and simulated velocity profiles summed over all outlet glaciers, see Supplementary Table 3); and all but one median correlation coefficient is larger than 0.74 (Supplementary Table 3). Observed horizontal surface velocities increase from near-zero at ice divides to B10 km per year in outlet glaciers, covering many orders of magnitude. This overall pattern is well represented in all calibration runs ( Supplementary Fig. 1) where flow arises due to the sum of longitudinal stretching and vertical shearing. For comparison, in a simulation where ice flow is only due to vertical shearing (see Methods), outlet glacier flow is basically absent ( Supplementary Fig. 2). This suggests that the spatial variability in flow can be explained to a large degree by the variability in ice thickness if membrane stresses are included.
With the calibrated model we performed simulations over a range of horizontal grid resolutions from 4,500 to 600 m. This allows us to address the impact of model resolution on simulation of flow in outlet glaciers, helping to set limits on the minimum resolution required to capture such flow.
Reproduction of flow pattern at the highest model resolution.
Of the 29 analysed outlet glaciers, the flow of 4 is characterized by low-surface slopes and low driving stresses ('ice-stream'-type), while all other glaciers flow through channels significantly deeper than the surrounding ice ('isbrae'-type) 13 . At our highest model resolution (600 m) the calibrated model reproduces the velocity structure of 'ice-stream' type glaciers slightly better on our cross profiles (median r ¼ 0.93) than 'isbrae' type glaciers (median r ¼ 0.88). Flow speeds in the transitional zone 100-300 km inland (speeds approximately 20-100 m per year) are generally underestimated (Fig. 1a,b).
Two glaciers, Ryder Gletscher (r ¼ 0.71) and Storstrømmen (ro0), are surge-type glaciers whose flow shows large variations over time due to hydrologically controlled sliding 14,15 , largely uncoupled from surface slope and ice thickness, which may have led to a reduced fit. In addition, Graulv and the Køge Bugt glaciers have correlation coefficients less than or equal to 0.65, substantially lower than the median (0.88). For these glaciers, Morlighem et al. 8 used kriging to derive ice thicknesses away from OIB measurements, while they used mass conserving interpolation techniques 7 for all other glaciers, enabled by better OIB measurement coverage. We hypothesize that the poor match between simulated and observed flow for these glaciers results from poorly constrained ice thickness knowledge. While glaciers are rapidly changing, our model is diagnostic and should reflect the flow for the current geometry. Results for other glaciers with large dynamic adjustments, such as Jakobshavn Isbrae, are very encouraging.
Role of model resolution. The impact of poor model resolution on the simulation of outlet glacier flow is illustrated by resampling the MO2014 subglacial topography and ice thickness to a range of model grid spacings (Fig. 2). To quantify the impact of improved model resolution, we performed a linear regression analysis of the r.m.s. velocity difference for each glacier. A simulation at the lowest model resolution (4,500 m) has an r.m.s. velocity difference 42% higher than at 600 m (Fig. 3c). Overall we find that r.m.s. velocity differences decrease with increasing resolution (goodness-of-fit r 2 ¼ 0.95). Sixteen glaciers have statistically significant (Po0.05) improving trends. All but one of these (Humboldt Gletscher) are fast-flow marineterminating glaciers (mean velocity 4200 m yr À 1 ) of the 'isbrae'-type. Jakobshavn Isbrae, Greenland's most prolific supplier of ice to the ocean, flows through a o10-km wide subglacial valley grounded well below sea level. A minimum resolution of B2 km is needed to resolve the valley geometry sufficiently to reproduce the high (42 km yr À 1 ) surface velocities (Fig. 2b, Supplementary Fig. 3). 'Ice-stream'-type glaciers and surge-type glaciers (Ryder Gletscher and Storstrømmen) and 'isbrae'-type glaciers with subglacial topographies derived by kriging (Graulv, Køge Bugt N and S; Supplementary Figs 3-5) do not benefit from increased resolution.
Role of ice thickness. To assess how much OIB measurements and mass-conserving interpolation methods have improved model simulations, we performed two simulations with an older ice thickness/subglacial topography data set (BA2001) 16 . The BA2001 data set is posted at 5,000 m, which we resampled to 600 and 4,500 m spacing. For these experiments, the total r.m.s. difference between modelled and observed flow of all 29 outlet glaciers are 49% (600 m) and 57% (4,500 m) higher than with the MO2014 ice thickness data set at 600-m resolution (Fig. 3, Supplementary  Table 4). At the same time, the reproduction of the flow structure is substantially better when using the MO2014 data set at a model resolution of 600 m (median correlation coefficient e r ¼ 0:58) than using the BA2001 data set (e r ¼ 0:58 and e r ¼ 0:25 at 600 and 4,500 m model resolution, respectively). High model resolution alone clearly is insufficient to reproduce outlet glacier flow (Supplementary Figs 6-8) Furthermore resampling data sets to resolutions higher than the nominal resolution (for example, resampling the BA2001 data set to 600 m) does not generate physical information that is not already captured by the coarse data. We note that, in some cases, such resampling may produce an apparent improvement. However, this improvement is an artifact of the analysis method (Supplementary Note 1).
Compared to 'isbrae'-type glaciers, 'ice-stream'-type glaciers are less sensitive to both model and data set resolution. Furthermore, the total r.m.s. velocity differences of 'ice-stream'-type glaciers is almost an order of magnitude smaller than for 'isbrae'-type glaciers (Supplementary Table 4).
The structures of three 'ice-stream' type glaciers that are well represented at the highest resolution (79North, Petermann Gletscher and Zachariae Isstrøm) are well represented at all grid resolutions. The improving trend of 'ice-stream'-type glaciers is not statistically significant (Fig. 2a). This suggests that all grid resolutions sufficiently capture lateral variations in ice thickness and surface slope, which are typically small in 'ice-stream' type glaciers. For 'ice-stream'-type glaciers accurate knowledge of ice thickness is less crucial to reproduce the flow structure than it is for 'isbrae'-type glaciers. Nonetheless subglacial topography still needs to be well constrained to properly treat advance and retreat scenarios.
Humboldt Gletscher. Humboldt Gletscher stands out in all aspects. Located in the northwestern corner of Greenland, the B90-km wide 'ice-stream' type glacier has a mean velocity o200 m yr À 1 . Humboldt's simulated flow structure agrees poorly with observations at all resolutions, shows the smallest improving trend, and has the lowest goodness-of-fit (r 2 ¼ 0.63). Simulated basal temperatures for Humboldt are very close to the pressure-adjusted melting point, and depending on grid resolution, basal motion may be activated or not in our simulations. The flow of Humboldt Gletscher is thus very sensitive to grid resolution and model parameters, unlike other glaciers studied. To explain this sensitivity, a more focused study might be warranted, possibly using time-dependent data assimilation and inverse methods such as in refs 17,18.

Discussion
The simulations presented in this study portray Greenland's flow structure in unprecedented detail; most large outlet glaciers considered are well-captured at grid resolutions o1 km. Our results demonstrate that spatial variability in flow can be explained in large part by the spatial variability in ice thickness. We find that inversion of surface properties for individual glaciers is not essential to reproduce the overall flow pattern. Using simple parametrizations of basal motion and subglacial hydrology, we find good agreement between simulated and observed spatial flow patterns while the reproduction of the velocity magnitude requires further improvements, especially in the transitional zone 100-300 km inland (speeds B20-100 mper year). Disagreement between observed and simulated speeds may arise from inadequacies in parametrizing basal motion and subglacial hydrology. Also not all heat sources that affect the viscosity of ice are accounted for in our model. For example refreezing of surface meltwater ('cryo-hydrologic warming') can soften the ice and enhance flow, which may be relevant in the transitional zone [19][20][21][22] . Furthermore the ice flux from the interior may be misrepresented because interior ice thickness remains less well characterized, but might be just as important to capture outlet glacier flow downstream. Finally, we calibrated only three parameters and only at the locations of our profiles, which may partly explain the reduced fit in the interior. To further reduce the gap between observations and simulations, future work must elucidate the physics behind the remaining misfit between modelled and measured flow across the spectrum of outlet systems, in the transitional zone and the interior. A quantitative, physically based understanding of the mismatch can then guide future modelling efforts.
Here we have presented diagnostic simulations; future work should address whether the observed temporal variability can be explained by our simple parametrizations or whether more physically based models of subglacial hydrology are needed. Basal motion can be highly variable locally on short (o1 year) time scales [23][24][25][26] , arising from variations in the hydrological system. Development of adequate models of subglacial hydrology, applicable on an ice sheet wide scale, is still at an early stage (for example, ref. 27). On a regional scale, recent efforts show promise in capturing the seasonal evolution in flow speeds resulting from changes in basal resistance caused by surface meltwater delivery 28 .
Here we have used surface velocities as our metric of success. It could be argued that ice discharge is a more suitable metric of success. However, ice thickness is subject to larger observational uncertainties than ice velocity, making discharge (the product of ice thickness and vertically averaged horizontal velocities) a weaker metric. We note that our results would look even more favourable if we used flux as a metric, with 20 glaciers showing correlation coefficients exceeding 0.85 ( Supplementary Fig. 9), and the fluxes of 6 marine-terminating glaciers being within observational uncertainty.
Clearly, large-scale observational efforts such as OIB improve our ability to model the parts of the system that are responsible for the bulk of the mass flux to the oceans. The simple metrics we introduced here help to quantify the impact of these additional measurements on whole ice-sheet simulations. Ice thickness measurements by many groups are ongoing and OIB, currently scheduled to continue until 2019, is expected to add many line kilometres of measurements by that time. Improvements such as those discussed here for Greenland would also be expected from additional investment in mapping subglacial topography for the Antarctic ice sheet. Our simulations have implications for efforts targeted at projecting twenty first century sea level rise. IPCC's Fourth Assessment Report called out the need to resolve the full stress configuration in ice-sheet models to simulate changes in outlet glaciers. Since then, ice sheet models have seen substantial improvements in their representation of flow physics (for example, refs 29,30). We find that models that resolve both membrane and vertical stress gradients are capable of reproducing the observed flow structure with high fidelity. In regions with large transverse velocity gradients, such as shear margins, the mismatch between observed and simulated flow may be further reduced by resolving additional components of the stress balance.
The large improvements we have attained using only spatially uniform parameters suggest that there is much to be learned about outlet glacier flow from comparing suites of glaciers to each other. By modelling each outlet glacier with the same spatially uniform parameters, modelled flow variations are dependent locally only on the shape of the bed and surface. Discordance with observed velocity, if not due to errors in these shapes, then may reflect departure from the conditions assumed by these uniform parameters. If groups of glaciers that experience similar surface melt and accumulation patterns show similar departures, this would indicate a common response that is not captured by present physics. In this way, it should be possible to limit the degrees of freedom in the model as the parameter space is enhanced to better capture spatial patterns of flow. This type of analysis would not be well bounded or even tractable if locally variable parameters had been used to match the observed flow patterns of each glacier.
Successful prognostic models will have to demonstrate that both spatial and temporal variability in flow can be replicated 2 ; here, we have shown significant progress in the former.

Methods
Outlet glaciers. Along cross-profiles of 29 outlet glaciers (Supplementary Table 2) we sample observed and modelled surface velocities every 250 m using bilinear interpolation, and all metrics are then calculated from those profiles. Profiles are drawn roughly perpendicular to observed flow directions, following flightlines whenever possible. Observed surface velocities 12 U s represent fall velocities of 2008; we assume that annually-averaged velocities U are 10 ± 5% higher than fall speeds, The flow component normal to the gate is given by U > ¼ n Á U, where n is the unit normal vector at each gate point such that n Á U is positive. We calculate the uncertainty in surface velocities, s, by adding 0.05U > to the 1-sigma error of the data set, U e : Ice sheet model. Simulations are performed with the open-source Parallel Ice Sheet Model (PISM), which is thermomechanically coupled, polythermal, and includes a hybrid stress balance model 10,31 . The hybrid scheme combines the Shallow Ice Approximation (SIA) 32 for vertical deformation and the Shallow Shelf Approximation (SSA) 33 for longitudinal stretching. The effective viscosity of glacier ice, Z, is given by where E is the flow enhancement factor, t e is the effective stress, A is the enthalpydependent rate factor (softness), and n is the exponent of the power law. The small constant e (units of stress) regularizes the flow law at low effective stress, avoiding the problem of infinite viscosity at zero deviatoric stress.
Boundary conditions. At the basal boundary, geothermal flux varies in space 34 , and a pseudoplastic power law 35 relates bed-parallel shear stress, s b , and the sliding velocity, u b : where t c is the yield stress, q is the pseudoplasticity exponent, and u 0 ¼ 100 m yr À 1 is a threshold speed. Here we use a simple model of basal hydrology that connects basal water pressure and the rheology of soft sediments ('undrained plastic bed') 36 . While this model was motivated by observations from Ice Stream B in Antartica 37 , there is growing evidence of soft sediments beneath the Greenland Ice Sheet [38][39][40][41] . Furthermore, theoretical work 42 suggests that hard beds, that is beds without soft sediments present, behave plastically at high basal water pressure. We assume that yield stress t c is proportional to effective pressure N ('Mohr-Coulomb criterion' 43 ), where f is the till friction angle and the effective pressure N is given by 37,27 : Here d ¼ 0.02 is a lower limit of the effective pressure, expressed as a fraction of overburden pressure, e 0 is the void ratio at a reference effective pressure N 0 , C c is the coefficient of compressibility of the sediment, W is the effective thickness of water and W max is the maximum amount of basal water. We use a non-conserving hydrology model that connects W to the basal melt rate _ B b (ref. 37): where r w is the density of water and C d ¼ 1 mm yr À 1 is a fixed drainage rate.
We prescribe the yield stress as a function of bed elevation 44 by prescribing the till friction angle f as a continuous function of the bed elevation, with f ¼ 5°for bed elevations lower than 700 m below sea-level, f ¼ 40°for bed locations higher than 700 m above sea-level, and changing linearly in between. While the bed elevation dependence is not derived strictly from first principles, outlet glaciers with bed elevations below sea-level have hydrologic systems tied to sea level at the terminus, and so run at high basal water pressure (that is, low effective pressure), and may be underlain by marine sediments 45 , which are weak 46 .
To assess the role of this bed elevation dependence, we performed a simulation where the yield stress is only a function of locally derived effective pressure, not dependent on bed elevation. In this simulation, fast flow is apparent in outlet glaciers; however, not of the right structure ( Supplementary Figs 10-12). The total (all 29 outlet glaciers) r.m.s. difference is 27% higher than the calibrated simulation, indicating that using an elevation-dependent yield stress provides a much improved match to the flow structure of most investigated glaciers.
Model initialization. The initialization procedure follows closely the 'flux-corrected paleo-climate' initialization method in ref. 47. During the last 5,000 years of the initialization we apply a flux-correction method to obtain an ice sheet geometry in closer agreement with measurements 47 . For computational efficiency grid refinements are made during the initialization procedure. The runs are started on a 18,000 m grid at À 125,000 years, then refined from 18,000 to 9,000 m at À 25,000 years, to 4,500 m at À 5,000 years, to 3,600 m at À 1,000 years, to 1,800 m at À 500 years, to 1,500 m at À 300 years, to 1,200 m at À 200 years, and finally to 900 m at À 100 years. The vertical model resolution is 20 m.
Model calibration. We calibrate three ice-sheet-wide scalar parameters, the flow enhancement factor E, the exponent of the sliding law q, and the exponent of the flow law for the SSA n, while all other parameters are PISM default. The calibration is performed at 1,500-m grid resolution.
To calculate the velocity field that corresponds to a given parameter combination, one could solve the non-linear system of equations of steady state using the observed geometry and standard direct or iterative solution techniques. Without smoothing applied to the observed geometry, such a velocity field would reflect any inconsistencies between the calculated thermodynamic and hydrological initial state and observations. Instead we attempt to follow the physical transient using a flux correction method 47 to iterate towards a geometry that is 'close' to the observed geometry and is consistent with the model and the applied forcing. Such a flux correction is used to minimize the difference between the modelled and the observed geometry. All experiments use the enthalpy field from the end of the initialization procedure at the corresponding grid resolution, except for 600-m grid resolution, which starts from the initialized state at 900-m grid resolution. All experiments are run for 100 years, sufficiently long to adjust to a given parameter combination. Flux correction does not respect observed ice thickness perfectly; especially in outlet glaciers near the grounding line, differences may locally exceed 100 m.
First we calibrate E by running the model in SIA-only mode and comparing simulated with observed surface speeds in slow-flowing (o20 m yr À 1 ) areas. The r.m.s. differences between modelled and observed surface speed range from 6 m yr À 1 (E ¼ 1.5) to 12 m yr À 1 (E ¼ 3; Supplementary Fig. 13). Except for E ¼ 3, all r.m.s. differences are smaller than the observational error of 9 m yr À 1 .
We choose E ¼ 1.25 to avoid overestimating the vertical shearing in slow-flowing areas.
Second, we calibrate the hybrid model by tuning the exponent of the pseudoplastic sliding law q, which connects bed-parallel basal shear stress and basal sliding, and the exponent flow law n for the SSA. For the SIA, we choose the default value n ¼ 3. To keep the parameter space manageable we test selected pairs of q and n (Supplementary Table 3). Here we use the r.m.s. difference along all profiles as our metric. The r.m.s. differences for all experiments are listed in Supplementary Table 3. By this metric, Experiment 7 with (q, n) ¼ (0.6, 3.25) scores highest. Therefore, all subsequent simulations are performed with (q, n) ¼ (0.6, 3.25).