Decadal-scale sensitivity of Northeast Greenland ice flow to errors in surface mass balance using ISSM

.


Introduction
[2] The Greenland Ice Sheet is a major source of future sea level rise, as it is vulnerable to warming trends, marginal ablation, and consequential retreat [e.g., Ridley et al., 2005;Mikolajewicz et al., 2007;Vizcaino et al., 2010;Huybrechts et al., 2011].While its response to warming is best understood on century and longer time scales, it is clear that the ice sheet is susceptible to dynamic ice flow variations that occur on much shorter (seasonal to decadal) time scales [e.g., Luckman et al., 2006;Rignot et al., 2008;Zwally et al., 2011;Bergmann et al., 2012].Several key processes have been linked to dynamic ice flow change, most notably penetration of surface runoff to the glacier bed [Zwally et al., 2002;Bartholomew et al., 2010] and grounding line retreat due to ice-ocean interactions [Rignot et al., 2010;Holland et al., 2008].Still, predictions of Greenland's response to shortterm climate change remain limited, due to an incomplete understanding of these processes, as well as the difficulty of accurately resolving surface mass balance (SMB) or the local surface accumulation minus runoff [Bales et al., 2009;Ettema et al., 2009;Burgess et al., 2010].Improved prediction of Greenland's ice flow requires a detailed understanding of SMB, basal hydrology, ice-ocean interaction, and the ice sheet's sensitivity to the uncertainties associated with each process.In this study, we specifically focus on determining Greenland's dynamic sensitivity to errors associated with SMB.
[3] Few modeling studies have investigated the influence of SMB on decadal-scale ice flow, in particular with respect to how it affects local mass flux.Most studies instead focus on predicting Greenland's long-term response to future warming scenarios, on the order of centuries to millennia, using coarse two-way model coupling between general circulation models (GCMs) and ice sheet models [e.g., Robinson et al., 2010;Vizcaino et al., 2010;Huybrechts et al., 2011].The ice sheet models typically used in these experiments are computationally efficient and practical for long-term coupling experiments, yet they often neglect longitudinal stresses.For modeling higher-velocity regimes, this assumption is disadvantageous, as it will break down where slip is dominant [Hindmarsh, 2004].Some studies, such as Koenig et al. [2011], have addressed the use of higher-order hybrid models that include longitudinal stresses to study Greenland's sensitivity to climate forcing.However, they focus exclusively on the ice sheet's evolution millions of years ago, and their climate forcing is spatially extrapolated from a regional-scale parameterization of accumulation and melt.Price et al. [2011] used a three-dimensional, higher-order ice flow model to explore Greenland's diffusive response to rapid changes in the marine termini of large outlet glaciers.The model does not account for changes to surface climate but does take advantage of higher-order physics to demonstrate how the sensitivity of fast-flowing outlets can be communicated to the ice sheet's interior.Most recently, Quiquet et al. [2012] conducted climatic steady state experiments using a hybrid Greenland Ice Sheet model.Their study focused on investigating Greenland's century to millennia-scale sensitivity to atmospheric forcing fields derived from various climate models.
[4] In this study, we investigate Greenland's sensitivity to decadal-scale climate variability.By quantifying the sensitivity of ice stream dynamics to errors in SMB, we aim to better understand the impacts of short-term climate trends on ice flow.Here we choose a high-resolution SMB forcing derived from Regional Climate Model (RCM) output.This product spans 171 years and is produced through a statistical reconstruction of RCM output that has been calibrated to match in situ observations and assessed for errors [Box et al., 2012;Box, 2013].Finely resolved spatial maps of surface conditions, like those derived from RCMs, are critical for accurately representing snowfall, especially in areas with steep topography or narrow channels [Glover, 1999;Ettema et al., 2009].For Greenland, only a few products of this sort are available due to the cumbersome process of interpolating sparse observations at meteorological stations onto an area of nearly 1.8 million km 2 [Ettema et al., 2009;Stocker et al., 2010;Hanna et al., 2011;Franco et al., 2012].Here we take advantage of such a product in order to investigate the sensitivity of ice flow to errors in SMB forcing.
[5] To determine how these errors propagate through an ice flow model, we use uncertainty quantification (UQ) methods to examine the calculated ice flux variability in a simulation of the Northeast Greenland Ice Stream (NEGIS) and its Zachariae Isstrøm outlet.The climatic sensitivity of this area is of great interest, as thinning has been observed at Zachariae Isstrøm since 1999 [Rignot et al., 2001], and an annual aerial survey based on Moderate Resolution Imaging Spectroradiometer (MODIS) images ranks this glacier third only to Humboldt and Petermann glaciers in terms of area of ice lost per year during the 21st century [Box and Decker, 2011].NEGIS ice flow is simulated on a high-resolution mesh by a two-dimensional (2-D) Shelfy-Stream Approximation (SSA) [MacAyeal, 1989] implemented within the Ice Sheet System Model (ISSM) [Larour et al., 2012c].The UQ methods provided by ISSM rely on Latin hypercube methods to derive uncertainties in model diagnostics, given specified errors in the model inputs.These capabilities have recently been used to quantify mass flux uncertainties in Pine Island Glacier, Antarctica, and to investigate how realistic errors in geothermal heat flux and ice thickness propagate through a static SSA model [Larour et al., 2012a].
[6] The SSA model is simple and less prohibitive than higher-order models, making it ideal for running samplingstyle UQ methods, like the analyses utilized here, in a realistic amount of time [Larour et al., 2012b].SSA is also based on equilibrium of longitudinal stresses, making it a viable choice for modeling fast-flowing regions such as the NEGIS, which has relatively flat bedrock with low driving stresses and a weak bed [Joughin et al., 2001].Past studies have used alternate methodologies to study the sensitivity of lower order Greenland models to surface parameters like surface accumulation [Heimbach and Bugnion, 2009] and summer temperature [Robinson et al., 2012].However, to our knowledge, no study has used similar methods to investigate how errors in SMB forcing propagate through a forward model that considers longitudinal stresses.
[7] In the first part of this study, we describe the ice flow model and the process used for UQ of SMB forcing.In the second part, we describe the model setup, including meshing and the data sets used for initialization and forcing.In the third and fourth parts, we present and discuss UQ results, focusing on how these results inform UQ of SMB forcing onto a forward model on decadal time scales.Finally, we conclude by discussing the dynamics captured by the model and reviewing what our results suggest about the ice flow and climate model improvements needed to achieve more accurate ice stream simulations.

Ice Flow Model
[8] ISSM is a thermo-mechanical ice flow model that relies upon the conservation laws of momentum, mass, and energy, combined with constitutive material laws and boundary conditions.The ISSM implementation of these laws and treatment of boundary conditions has been described in Larour et al. [2012c].ISSM relies on the finite element method, which allows for simulation of ice flow on an anisotropic mesh, a feature advantageous for modeling areas with strong acceleration or within narrow channels at finer resolutions.In terms of initialization, we rely on inversion methods to ensure that the ice flow model matches the present-day configuration [Morlighem et al., 2010;Larour et al., 2012c].Below, we review the ice sheet model components pertinent to this study, namely the mass transport and mechanical models.
[9] For this study, only SMB is forced through time.No thermal model is invoked here, so ice temperatures remain fixed through time and averaged vertically on a 2-D grid.The model does not included bedrock deformation or migrating boundaries; thus, the bed geometry and the ice sheet area remain fixed.Additionally, parameters for basal drag and all other boundary conditions are held constant throughout the simulation.This ensures that all ice responses are forced only by the spatial and temporal evolution of the SMB.
[10] Changes to SMB forcing are directly communicated to the ice flow through the mass transport model.This model is driven by mass conservation, which relates the rate of ice thickness change, the divergence of ice flow, surface mass balance, and basal mass balance as follows: where H is ice thickness, v = (u, v) is the depth-averaged horizontal velocity, P M s is the surface mass balance (m/yr ice equivalent), and P M b is the basal mass balance (m/yr ice equivalent).Note that basal hydrology is not simulated here; thus, the basal mass balance is set to zero.
[11] To calculate velocities, we model mechanical ice flow with a 2-D SSA [MacAyeal, 1989].We refer the reader to Larour et al. [2012c] for a full description of the implementation of these equations within ISSM.The SSA is based on the full-Stokes equations, assuming that vertical shear and bridging effects are negligible, reducing the momentum balance equation to a two-dimensional system of equations.
[12] At the ice-bedrock interface, we apply basal drag empirically after Paterson [1994]: where v is the horizontal velocity vector, N the effective pressure of the water at the glacier base, b the basal shear stress, and ˛the basal drag coefficient.The evaluation of effective pressure requires a full hydrological model, which is a work in progress [Larour et al., 2012c].As a first order approximation, we use N = g( H + w z b ) [Paterson, 1994], where and w are the density of ice and water respectively, g is gravity, and z b is bedrock elevation.The value of H, and therefore N, varies as the ice sheet geometry changes with time.At the ice-sea water interface, depth-integrated water pressure is imposed.On all other boundaries where stresses are not specified, single-point velocity constraints equal to the observed velocities [Rignot and Mouginot, 2012] are applied.

Uncertainty Quantification
[13] Our UQ methods are based on the Design Analysis Kit for Optimization and Terascale Applications (DAKOTA) software [Eldred et al., 2008], which is embedded in ISSM.These capabilities were previously established and validated by Larour et al. [2012aLarour et al. [ , 2012b]].ISSM offers two types of UQ methods: sampling analysis and sensitivity analysis.Sampling analysis is typically used to quantify the uncertainty propagated through a forward model in response to errors in model inputs.This analysis relies on repeated execution of the same model, where input variables are perturbed within specified distributions simultaneously at each partition for each individual run.Sensitivity analysis, on the other hand, quantifies how the location of errors impact mass flux through a specified flux gate.Sensitivities (change in output with respect to a change in input) are calculated by perturbing each partition one at time and assessing how changes within the specified partition affect model output.These analyses are carried out on equal-area partitions of the model domain.Partitioning is accomplished using the CHACO Software for Partitioning Graphs [Hendrickson and Leland, 1995].

Sampling Analysis
[14] For sampling analyses, we run the same forward simulation many times in order to assess the effect of random errors in SMB forcing.For each individual run, input variables are perturbed by different amounts at each partition.Input values are perturbed randomly within a prescribed range (described by a statistical distribution, e.g., normal or uniform), separately for each partition.For  et al., 2012;Box, 2013].Outlined in gray is the NEGIS regional domain (ISSM-NEGIS) of the 22 year run used in uncertainty quantification.This domain is extracted from the larger continental domain after completion of a 150 year relaxation.
example, a normal distribution for a particular partition is fully described by an average, , and a standard deviation, .By definition, normal distributions cluster around and decrease towards the tails, in a Gaussian bell curve ranging at ˙3 .On the other hand, the uniform distribution places greater emphasis on values closer to the tails, where probability of occurrence is equal for any given value within a radius r of the average, .Here the value of r is set equal to the value of 3 .Comparison between responses of both uniform and normal distributions offers advantageous insight into the model sensitivity to extreme forcing.Therefore, in this study, we sample with both normal and uniform distributions.
[15] At the beginning of a given simulation run, a random percentage perturbation P i , is chosen for each partition i.Its value falls within a distribution (-3 i Ä P i Ä 3 i ), where Every e i is defined separately for each individual partition i, according to the standard errors plotted in Figure 1.P i remains constant for every partition over the duration of each 22 year simulation.At every time step t, the forcing SMB value within a partition i is perturbed by P M s (t)P i .Resulting statistics are calculated for the output responses at the conclusion of the simulation.This includes means, standard deviations, and cumulative distribution functions (CDFs).Generation of the P i values relies on a binned Latin hypercube sampling (LHS) algorithm [Swiler and Wyss, 2004], which has the advantage of forcing samples into the tails.
It is also a more efficient method of sampling for a given level of statistical accuracy [Larour et al., 2012b].In this study, our input forcing is SMB, and the output responses of interest are mass flux through prescribed gates (Figure 2b), maximum velocity, and ice volume.

Sensitivity Analysis
[16] To perform a sensitivity analysis within our domain, ISSM carries out an assessment of the local sensitivity of the output responses to perturbations in the input forcing.
Using the established ISSM framework, we assess temporally integrated errors in model output as ice thickness and velocity evolve over time.These errors are affected by the feedback between the mass transport and mechanical models as the ice flow responds to forcing.To investigate these processes, we quantify how the location of errors impact mass flux through a specified flux gate by imposing a 0.1% change in SMB to a different partition per each individual run, resulting in one 22 year simulation per partition.Mass flux responses are determined at the completion of the simulation.We refer the reader to Larour et al. [2012b] for a discussion about the propagation of errors through a static ice flow model, associated assumptions, as well as the calculations of mass flux output and response sensitivities using a forward model.
[17] Sensitivity values within each domain partition, Â i , represent the magnitude of output response to input perturbation.In this study, SMB error inputs are scaled to a value of 1, so sensitivities are given in units of mass flux response.We additionally calculate the scaled sensitivities, SS i , which provide nondimensional sensitivity values within the n domain partitions, where: Scaled sensitivities nondimensionally represent the relative contributions of sensitivities throughout the domain.Sensitivity values over the entire domain sum to unity, so they can be used to compare the relative contribution of model inputs.

Forcing
[18] The SMB time series is interpolated from 5 km yearly grids covering the entire Greenland Ice Sheet [Box et al., 2012;Box, 2013].This product is a statistical reconstruction of surface accumulation and melt during hydrologic years (October through September) 1841-2010 and is driven by output from the Polar NCAR/Penn State Mesoscale Model (PMM5) [Box et al., 2012;Box, 2013].Its advantage is that it combines the spatial detail of a regional model with the absolute accuracy of in situ data, while spanning a much longer time period than available from other RCMs (which typically begin in 1958).
[19] The melt reconstruction is based upon regression between PMM5 temperature output and long-term surface air temperature records [Box et al., 2012;Box, 2013].PMM5 melt water production is computed from mean monthly near-surface air temperatures using a positive degree day scheme after Braithwaite [1985], Reeh [1989], andFausto et al. [2009].Separate degree day factors are used for snow and bare ice surfaces.A simple meltwater retention scheme yields residual runoff and bare ice area.In the ablation zone, the SMB reconstruction has a root mean square error (RMSE) of 0.6 m/yr water equivalent.Therefore, for the purposes of this study, standard error in the ablation zone is set to a value of 0.6.It is important to note that, in areas of extreme melt (ablation > 4 m/yr), the reconstruction overestimates ablation by about 2 m water equivalent, though in this study, we will not take these extremes into account.
[20] Accumulation reconstruction is based upon regression between PMM5 accumulation output and ice core measurements after Box et al. [2009].There are a total of 82 ice cores.For a given year, the available ice cores contribute separately to the reconstructed values, and the cores that correlate best with PMM5 output are given the strongest contribution factor.The new accumulation grids are improved over Burgess et al. [2010], as the interannual amplitude of accumulation has been matched to ice core measurements.Accumulation standard errors are equal to one standard deviation of regression residual.For this study, standard error in the accumulation zone is set equal to these standard deviations [Box, 2013].The combined standard error map is plotted in Figure 1, which has a distinct patchwork appearance.This feature is an artifact of the regressions, as values of single meteorological stations drive the reconstruction at each grid cell.
[21] SMB is imposed through a one-way coupling scheme [Huybrechts et al., 2004;Gregory and Huybrechts, 2006;Stone et al., 2010;Schlegel, 2011], which simplifies our model and is advantageous in deconstructing ice flow sensitivity to errors in SMB.A two-way coupling scheme can often be complex, as it requires downscaling techniques in order to translate between low-resolution climate model output (1°-3°) and high-resolution ice flow model input (1-40 km) [Bougamont et al., 2005;Vizcaino et al., 2010;Robinson et al., 2012;Helsen et al., 2012].One-way coupling allows us to use an RCM-derived product, which provides forcing at a spatial resolution on the same order as the ISSM mesh, while eliminating the need for complex coupling routines.

Inversion
[22] Initialization is carried out for the entire Greenland Ice Sheet, using an inversion technique [MacAyeal, 1993] following Morlighem et al. [2010], where we determine the basal drag coefficient required to match InSAR surface velocities from Rignot and Mouginot [2012] (Figure 2c).Prior to inversion, initial properties of the ice sheet (i.e., geometry, velocity, and surface boundary conditions) are prescribed according to observations.The Greenland Ice Sheet model anisotropic mesh is made up of 73,320 elements with a resolution ranging from 1 km at high-velocity outlets to 15 km at the ice divide.The bedrock geometry is initialized using 5 km gridded bedrock data from Bamber et al. [2001] and gridded data for the Jakobshavn, Petermann, Helheim, and Kangerdlugssuaq Glacier regions derived from NASA's Operation Ice Bridge measurements [Leuschen et al., 2011], available from the Center for Remote Sensing of Ice Sheets (CReSIS) [Gogineni, 2012].The ice surface is after Scambos and Haran [2002], and SMB forcing is provided by Box [2013] (Figure 2d).Prior to inversion, ice temperatures are inferred by solving for thermal steady state using present-day boundary conditions including surface temperatures from Ettema et al. [2009] and geothermal heat flux from Shapiro and Ritzwoller [2004].The thermal model is as described in [Larour et al., 2012c] and includes advection and diffusion in all three directions.3-D transport velocities are calculated using a stress balance model, forced by observed surface velocities [Rignot and Mouginot, 2012] and an initial estimate of basal friction.Ice rheology is calculated by vertically averaging the thermal steady state ice temperatures.During the inversion, the prescribed thermal regime remains constant.All ice is considered grounded, and the ice front position is imposed and fixed in time.

Relaxation
[23] The ice sheet resulting from inversion is not in steady state.Evolution of the ice sheet from this initial configuration leads to changes in ice volume as ice sheet thickness and velocity continue to adjust towards a state where total ice sheet mass balance is zero.To limit this effect, we allow the model to adjust to steady state before attempting to run forward through time.To do so, the model is forced with a constant SMB until thicknesses and velocities are no longer changing over time (i.e., discharge is equal to total SMB forcing).It takes 6000 years to reach this steady state configuration.As forcing, we use average SMB from the years 1971-1988, a period in which Greenland's total mass balance was close to zero [Rignot et al., 2008].Ice temperature fields are held constant throughout the simulation.We choose this simple procedure instead of running the model through several glacial cycles.It would be inappropriate to run the model through glacial cycles since we spin-up using observed velocities and present-day ice sheet conditions.In addition, ISSM does not allow for ice sheet boundary migration, so ice sheet advance or retreat would not be accurately simulated during a paleoclimate spin-up.Therefore, we adopt a simple spin-up procedure that isolates responses to perturbed forcing, as we are strictly interested in the sensitivity of the ice sheet to changes in SMB.After relaxation, the new ice sheet geometry differs most dramatically from the observed initial condition along the margins.In these areas, the relaxed ice sheet is on average 200 m thicker and 200 m/yr slower than presently observed.These changes are responsible for a 3.7% total increase in volume over 6000 years.
[24] After achieving steady state, we run the resulting ice sheet model through the years 1841-1988 using Box's yearly SMB time series [Box, 2013].Ice temperature fields are held constant throughout the simulation, ensuring that all responses are forced only by changes to SMB.The model state at year 1988, dictated by ice thickness, bedrock elevation, and velocities (Figure 2b), becomes the basis for our UQ.We acknowledge that this spin-up procedure is idealized, and as a result, we do not expect the relaxed ice sheet to have the same geometry as the actual ice sheet in the year 1988.Likewise, we do not expect the simulation runs to reproduce specific events in the past or predict future events.This is because of the following: (1) the actual presentday ice sheet is the result of past glacial cycles, responsible for large fluctuations in temperature and precipitation over hundreds of thousands of years; (2) our configuration is dependent on a number of independent observational data sets that are not necessarily compatible with each other [i.e., Seroussi et al., 2011;Morlighem et al., 2011]; and (3) ISSM does not simulate ice sheet hydrology or ice-ocean interactions.However, since our goal is to use an ISSM configuration that resembles the present-day Greenland Ice Sheet in order to determine the relationship between ice flux and SMB, the spin-up procedure implemented is reasonable.
[25] Before running the UQ simulations, we extract from the large continental domain a smaller regional domain, seen  2b, maximum velocity (Vm), and ice volume (Vol).We sample a 22 year run 200 times.All runs are forced with the same base SMB time series, but for each sample the forcing is perturbed spatially by random values that fall within a range of ˙3e (Figure 1).Red and blue correspond respectively to uniform and normal uncertainty distribution samplings.Output means ( ) and distribution range () relative to (in %) are provided.
outlined in Figure 1.This new domain, hereafter referred to as ISSM-NEGIS, which encompasses the NEGIS drainage basin, has a mesh of 3711 elements (Figure 2a).At the iceice boundaries of ISSM-NEGIS, ice velocities are forced with single-point constraints set equal to observed InSAR surface velocities [Rignot and Mouginot, 2012], while velocities within the domain and at the ice-water boundaries are allowed to evolve.The above extraction method assumes that ISSM-NEGIS captures the entire NEGIS drainage basin and that the velocities at the divide will be minimally affected by changes to SMB within the domain.This is a valid assumption over the short decadal-scale time period used in uncertainty quantification.

Uncertainty Quantification Methods
[26] After extraction, we force the smaller domain with the remaining 22 years of SMB forcing .During these last 22 years, we perform sampling and sensitivity analyses.The sampling analysis standard deviations are determined by the SMB standard errors (Figure 1).
[27] Eight flux gates are positioned throughout the ISSM-NEGIS domain.Mass flux is calculated within the model a Bold entries are significant above the 99.5% confidence interval.Gate 8 mass flux correlates negatively with gate 7 mass flux.Ice flux responses at gates 1 and 5 correlate positively to each other, and they both correlate negatively with the response at gate 6.
as described in Larour et al. [2012b].Mass fluxes through the gates are computed simultaneously throughout the model run, and here we use them as a metric to determine how errors in SMB alter the mass being advected through tributaries of the glacier's outlet.The model also keeps track of the maximum velocity and the ice volume of the NEGIS model domain.
[28] The model time step is 1.825 days, or 200 time steps per year.Model runs are launched on the NASA Advanced Supercomputing (NAS) Pleiades cluster, on 60 cpus, for approximately 80 h.For the sampling analysis, the domain is partitioned into 1000 equal-area domains, and standard error is averaged within each partition.Results reflect the distribution of output from 200 independent samples of SMB.For the sensitivity analysis, the basin domain is split into 500 equal-area partitions, and each partition is assessed for influence upon the eight mass flux gates (Figure 2b) given a 0.1% perturbation in SMB (Figure 2d).

Sampling Analysis
[29] Sampling analysis allows us to quantify the uncertainty propagated through a forward model in response to random errors in model inputs, in this case SMB forcing over a 22 year period.In Figure 3, we present the statistical distributions of sampling results for all eight flux gates as well as the domain volume and maximum velocity.In Table 1, we list the correlation coefficient between the sampling responses of all gates.
[30] Results are presented from two types of sampling statistical distributions: normal and uniform (Figure 3).Both distributions yield similar results, suggesting that the model tends to smooth out ice velocities, even when forced with extreme values of SMB.Gate 8 is the exception.Its uniform distribution is flatter, indicating a much higher sensitivity to extremes.This is also evidenced by the distribution's distinct double-humped shape, which indicates that mass flux is responsive to small perturbations from the mean.The double-hump is most dramatic for gates 4 and 8, but it is clear that all distributions share this characteristic.Further analysis reveals that changes to surface slope are responsible for this specific shape.Indeed, samples promoting steeper slopes cluster to the right of the mean, and samples promoting flatter slopes cluster to the left.Gates 4 and 8, located in the steepest area of the ice stream, are most sensitive to this effect.
[31] Furthermore, all the sampling distributions, with the exception of volume (skew > 2), have little skew (skew ranges from -0.35 to 0.27) and small tails, emphasizing the smoothing effect.The volume distribution, in contrast, has a long positive tail, which suggests bias towards mass accumulation.A slight skew is somewhat expected, as errors are sampled over the entire domain, and the area of accumulation (positive SMB) is much larger than the area of ablation (negative SMB) (Figure 2d).Much of the ice in the accumulation area flows slowly, so transport out of the system is delayed, leading to accumulation storage in higher elevations.This feature is further enhanced by the SMB forcing time series itself, as the distribution of yearly SMB is positively skewed and P i is held constant for all partitions throughout the simulation.However, we find that the most   extreme values in the positive tail are caused by dynamic changes to ice flow, forced by large errors in the most southern branch of the NEGIS (terminating in Storstrømmen glacier, 77°N, 22°W).This error is defined within a partition that coincides with a broad equilibrium line area (Figure 2d, blue shading).Here SMB values are small and negative, but standard errors are large.When a large and negative P i is sampled, ice velocities slow dramatically as the slope is flattened.Storstrømmen is effectively shut down, resulting in massive ice storage and consequential volume growth.The opposite forcing (a large positive P i ) does not result in such an extreme result.Instead, a dramatic increase in velocity thins the ice, preventing the ice stream from sustaining high mass outflux.
[32] Generally, we find that average mass fluxes are comparable to observations even after relaxation and spin-up [Rignot et al., 2008], ranging from 0.61 at the smallest tributary to 13.43 Gt/yr through the main upstream trunk of the NEGIS.The mass flux distribution uncertainties range from 1.6% to 16.5% for both normal and uniform sampling distributions.The distributions for maximum velocity are similar to those presented for gate 8 (Figure 3, Mf8).This is not surprising considering that maximum velocity occurs in close proximity to the mouth of the Zachariae Isstrøm outlet glacier.In the narrow outlet, partition resolution is coarse compared to the mesh resolution, so we would expect mass flux downstream of gate 8 to respond uniformly.The maximum velocity uncertainty is over 11% of its mean value.This uncertainty is significant, as it exceeds the observational uncertainty of InSAR velocities (about 1% in the high-velocity areas of the NEGIS domain) [Rignot and Mouginot, 2012].Volume uncertainty is 0.21% and 0.38% of total mean volume, equivalent to 1.7 and 3.0 mm of sea level  equivalent over the sampling period.As discussed above, the volume distribution is skewed due to a dynamic shutdown of the southern branch of the ice stream.Thus, it is important to note that while percent uncertainty (relative to total volume) is small, mass flux in this area can be significantly perturbed due to errors in SMB.
[33] Results in Table 1 offer evidence that ice flow within a 60 km radius is interconnected.For instance, the mass flux at gate 8 (the gate closest to the NEGIS outlet) is negatively correlated with the mass flux of gate 7 (which lies about 30 km upstream).The area between gates 7 and 8 lies in a regime of high-velocity flow and corresponds with the extent of the interannual migration of the equilibrium line (zero SMB, blue shading in Figure 2d).Gate 8, which lies within the ablation zone, also lies within an area of highvelocity flow approaching the glacial margin.Mass flux here is strongly correlated to the domain's maximum velocity (r = 0.75) as well as the SMB immediately upstream from the gate (r = 0.79), offering evidence that coupling exists between the ice flow at the glacial margin and ice flow tens of kilometers upstream.We find similar interconnections between mass flux gates further upstream, as ice flux through gates 1, 5, and 6 are significantly correlated (though less than gates 7 and 8).These gates lie within 60 km of each other.Ice flux through gates 1 and 5 correlate positively, and they both correlate negatively with ice flux through gate 6.The correlation between gate 1 mass flux and the mass flux in the ice stream's main trunk is noteworthy.This result suggests that ice flow in a tributary may significantly influence the ice dynamics of an ice stream.

Sensitivity Analysis
[34] Sensitivity analysis results offer a spatial perspective of how small perturbations to SMB can propagate on a decadal time scale.Sensitivities themselves are relative, so we can compare them over gates and upstream/downstream of the same gate.In Figures 4, 5, 6, and 7, we present the sensitivity analysis results for mass flux at gates 1-8 (Figure 2b).In Figures 4 and 5, we plot the sensitivities of mass flux to perturbations of 0.1% of the SMB forcing (Figure 2d).In Figures 6 and 7, we plot the corresponding scaled sensitivities.For each flux gate, Table 2 summarizes the radius of influence for different orders of magnitude of scaled sensitivities.The radius of influence is the maximum distance between each flux gate and all partitions that have scaled sensitivities less than a specified value.
[35] Results of sensitivity analysis clearly illustrate that positive upstream perturbations and negative downstream perturbations in SMB forcing promote increases in mass flux (Figures 4 and 5).These results suggest that mass flux is strongly tied to local surface slope.Gate 8 is an exception to this observation because of its proximity to the glacial margin, where the ice interfaces with water at sea level.
Here small positive perturbations in SMB upstream and downstream from gate 8 both promote an overall steepening downstream, positively influencing mass flux.Gates 4 and 7, located just upstream of gate 8, are also located on steep terrain.Their sensitivities show a strong preference for flow downstream, as thickening promotes a much weaker a Values are given for multiple orders of magnitude of scaled sensitivities (log 10, ranging from 10 -1 to 10 -10 ).Small influences can be found as far away as 200 km from the flux gate.On average, SMB has a first-order influence within a radius of 40 km.
absolute response to thickening downstream (faint blue shading) than upstream (red shading).These results further illustrate that the behavior of a gate is dependent on surface slope, suggesting that steeper terrain is more resistant to a decrease in mass flux due to SMB-forced flattening.In agreement with results presented in Table 2, sensitivity plots also suggest a negative correlation between the mass flux through gates 7 and 8.We find that the area immediately upstream from gate 8 influences mass flux through both gates, positively influencing flow through gate 8 and negatively influencing flow through gate 7 (Figures 5 and  7, gates 7 and 8).Another noteworthy feature is the general shape of the sensitivity contours.Close to the gates, there is little lateral influence, while sensitivities expand laterally with distance from the gates.These results highlight a decoupling of lateral influence with proximity to the gates, indicating a strong dependence on the direction of ice flow both upstream and downstream.
[36] By mapping scaled sensitivities, we gain a better understanding of the spatial extent of SMB forcing influence of mass flux.Scaled sensitivities decrease exponentially with distance to the flux gate, decreasing by 3 orders of magnitude at a distance of about 100 km (Figures 6 and 7).In fact, small influences can be found as far away as 200 km from the flux gate (Table 2).These are significant distances, suggesting that errors far upstream and downstream of a gate could affect local mass flux.In addition, results presented in Table 2 indicate that SMB has an average first-order influence within a radius of 40 km.In agreement with sampling analysis, scaled sensitivities offer further evidence that mass flux from an upstream tributary depends upon the flow of the ice stream's main trunk.In Figure 6, gate 1, it is clear that mass flux through gate 1 is influenced by ice flow downstream.The downstream contours turn along with the direction of ice flow, suggesting that influence is transferred along flowlines.Repetition of these sensitivity experiments using a larger SMB perturbation (10%) yielded similar conclusions (not shown), suggesting that these results are indeed robust and that errors do propagate linearly through the model.

Discussion
[37] Sampling results indicate that the mass flux distributions are not skewed and that the model velocities respond smoothly to noisy SMB input.These results suggest that the simulated mass flux is robust, especially for gates upstream of the main glacier outlet (gates 1 through 7).This observed linear response of ice flow to perturbations in SMB is similar to results presented by Larour et al. [2012b].Our results support their conclusion that the SSA model has a tendency to dampen large perturbations placed upon the system.This study extends their findings to a forward ice flow model.Results suggest that the mechanical model dampens the model response and serves as a filter even when run forward through time.When we compare the results from the normal with the uniform distribution, we notice that the results are similar for gates 1 through 7, further illustrating the filtering effect of the SSA model.Gate 8 is a unique exception, as it lies near the equilibrium line altitude (ELA) and within the region of maximum ice velocity.In response to uniform sampling, it exhibits a double-humped mass flux distribution, and its tails are not reduced like the other gates.These results illustrate that gate 8 is sensitive to perturbations in local SMB and that it is responsive to extremes, suggesting that narrow, fast-flowing outlets are more vulnerable to changes in SMB.Overall, we find that simulated mass fluxes over the entire NEGIS are robust and well behaved, though fairly sensitive to the current errors in SMB, especially near the margins.Thus, current errors in SMB forcing may not be small enough for simulating the ice flow of large ice streams in Greenland or Antarctica.Mass flux may be well behaved due to the lack of positive feedback associated with melt water hydrology and ice/ocean calving responses not modeled in this study.It is believed that these dynamic processes play an important role in the recent acceleration of Greenland's large outlet glaciers.Inclusion of these processes is necessary for improved model predictability of ice sheet response to rapid decreases in SMB.In the future, exploitation of additional ISSM capabilities (e.g., 3-D, full-Stokes equations, grounding line migration, integration of a more advanced computational sampling solver) will determine if more advanced physics would yield significantly different discharge sensitivities to SMB input.
[38] Overall, we find that changes to flow are not locally isolated.In fact, every tributary plays a part in defining ice stream flow, even those well upstream from the main outlet glacier.Gate 1 scaled sensitivities most strikingly illustrate this effect (Figure 5).Mass flux through gate 1, which spans a major tributary, is clearly affected by flow in the main ice stream.Significant correlations between ice fluxes at gates 1, 5, and 6 (Table 1) further emphasize the linkage between local ice flow and nonlocal ice flow 60 km away.Sensitivity results also reveal that even though the areas of influence above and below the flux gates are equal in size, their shapes are dependent on speed and direction of flow upstream and downstream of the gates.These results offer evidence of bias towards positive sensitivities, especially in proximity to the gates and in the direction of ice flow.We find that meandering ice flow clearly influences the sign of the sensitivity values, as the border between positive and negative sensitivities is not restricted to the gates (e.g., Figure 4, gate 3 and Figure 5, gate 5).In areas where there is a shift in the direction of ice flow or a convergence of different branches of the ice stream, a downstream increase in slope is the most preferential condition for increased mass flux.In order to capture these details and accurately simulate meandering ice flow, increased resolution in outlets is critical.Local mass flux has a complex sensitivity to nonlocal ice flow, suggesting that flowline models, which rely only on perturbations of ice thickness upstream and downstream of the gate, would not adequately represent a full response to a shift in climate trends.
[39] Larour et al. [2012b] suggest that ice thickness is a highly influential parameter in the simulation of ice stream mass flux using a mechanical model, largely due to the influence of thickness on local driving stress.Driving stress ( d ) is defined as follows: where ˛s is the local surface slope.Here thickness and surface slope together exert a direct influence on ice velocities through this driving stress relationship and the mechanical model (SSA).Indeed, SMB communicates with the simulated ice stream mass flux primarily through the mechanical model, where promotion of a steeper surface induces faster flow.Likewise, promotion of a flatter surface slows the ice flow.This is especially true within the ISSM-NEGIS, where bedrock is relatively flat with low driving stress and a weak bed [Joughin et al., 2001].As evidenced by the sensitivity results, mass flux does strongly depend on driving stress.In fact, the latter dominates ice flux response to perturbations in SMB forcing, as upstream increases and downstream decreases in SMB correspond with greater mass flux (Figures 4 and 5).For gates 1 through 7, it is clear that adjustments to driving stress can impact mass flux, indicating ice flow responds directly to changes in surface slope.Sampling analyses reveal that surface slope steepening does promote an increase in mass flux as evidenced by the distinct double-humped shapes of resulting mass flux distributions (Figure 3).Because our simulation depends on coupling between a mechanical model and a mass transport model, we may expect UQ analysis to offer insight into additional, secondary dynamic responses.For example, by equation ( 1), slope-driven changes to velocity affect mass transport through feedback to the flow divergence term (r Hv).Though this relationship suggests that SMB affects mass flux in a highly nonlinear way, we find no evidence of a secondary response or additional feedback in our results.Instead, support the assessment that, over decadal time scales, driving stress responses dominate ice flow adjustment to trends in SMB.Indeed, we find that non-uniform perturbations to SMB are capable of altering the surface slope and the driving stress, resulting in changes to local ice flow.
[40] Furthermore, results offer strong evidence of tight coupling between the ice flow at the glacial margin and ice flow at least 50 km upstream.Though sensitivity results suggest that small positive perturbations in SMB both upstream and downstream from gate 8 promote an increase in mass flux, the area directly upstream is of particular interest.This is specifically because the sign of SMB forcing here dictates the equilibrium line altitude (ELA) (Figure 2d, blue shading).By definition, SMB values in proximity to the equilibrium line are small, resulting in large percent uncertainties.This suggests that errors here play a significant role in dictating where the ice sheet will lose mass and where it will gain mass.Since the equilibrium line also happens to correspond with an area of high-velocity flow, its position is highly influential on both domain maximum velocity and the domain's total mass flux (gate 8).This is also demonstrated by the extreme positive skew in the volume distribution.Here large changes to surface slope near the ELA are responsible for shutting down the southern branch of the ice stream, resulting in mass storage and large volume increase.In fact, scaled sensitivities offer evidence that the ELA's influence on ice flow extends tens of kilometers upstream from the glacial margin.We observe that the ELA can shift upstream when i > 1 in the upstream partition and a large and negative P i is generated.It is evident that, even though this upstream shift in equilibrium line tends to decrease flux downstream, it enhances the driving stress tens of kilometers upstream, resulting in faster flow upstream.Since the downstream response (e.g., Figure 3, Mf8) is 4 times larger than the upstream response (e.g., Figure 3, Mf7), it is the decrease in outflux that dominates during the short time scale sampled in this study.
[41] We find that gates within tens of kilometers from each other are interconnected, as they share an area of high influence that sits in between (Figures 6 and 7).In proximity to Zachariae Isstrøm, this area of influence establishes communication between the outlet glacier and upstream continental ice, as evidenced by the significant correlation in ice flux between gates 7 and 8 (Table 1).These results illustrate a mechanism for upstream propagation of marginal thinning or draw down [Payne et al., 2004;Howat et al., 2007;Price et al., 2011], a process that is hypothesized to influence regions hundreds of kilometers upstream over periods longer than the 22 years explored in this study.In our simulation, we find that draw down is slowed by the negative feedback between SMB and the domain's total mass flux (gate 8).Indeed, we do not observe a marginal speedup in response to local increases in marginal ablation, though observations suggest that draw down occurs within decades in steeper outlet glaciers of the Greenland Ice Sheet [Howat et al., 2007].This is most likely due to the fact that processes known to enhance dynamic ice response (e.g., hydrology, grounding line migration, ocean dynamics, and interannual forcing) are not simulated here.Still, we observe that ice flow upstream is interconnected with ice flow hundreds of kilometers away.These results offer evidence that flow of the interior ice is coupled to ice flow at the ocean interface, through ice dynamics, suggesting that a mechanism for positive feedback response over longer time scales not only exists within the modeled system but is sensitive to perturbations in SMB.
[42] On average, influence on mass flux can extend beyond 200 km upstream and downstream of a gate (Table 2).Local changes within 40 km both upstream and downstream are found to be the most influential, with importance decreasing exponentially with distance from the gate.Further investigation (not shown) reveals that this is indeed due to the dynamics of the mechanical model, as sensitivity experiments forced with constant velocities suggest that the radius of influence for advection due to mass transport is less than 40 km.Thus, in the case of our coupled model, the mechanical and the mass transport model are together responsible for first-order uncertainties within this radius, while beyond 40 km, the mechanical model dominates.As the mechanical model is capable of propagating SMB errors hundreds of kilometers upstream and downstream, it is clear that accurate SMB is not just important near the mass flux gates.We find that accuracy throughout the entire domain is critical for realistic simulation of branching ice flow.Additional detail is especially important near the ELA where SMB values tend to be small.Our results illustrate that errors in SMB, especially in proximity to the equilibrium line, can alter mass flux estimates by over 16%.Though current error margins are relatively acceptable for large ice streams, these results indicate that attention must be placed on increased resolution and accuracy of SMB forcing.
[43] As this study reveals, first-order sensitivity of mass flux to errors in SMB falls within a radius of 40 km (Table 1), demonstrating that accurate ice flow simulation in this region requires SMB forcing to be provided at a spatial resolution of at least 40 km.Extending these results to the rest of the ice sheet is not straightforward, since the NEGIS region is unique.This large area (700 km in length) of wide and flat, fast-flowing ice has relatively low accumulation and little melt.In the Northern Greenland, SMB gradients are similar and slopes are flat; thus, we expect sensitivity experiments would yield similar results in other fast-flowing outlets in this area.In the southern half of the ice sheet, however, ablation is much stronger, leading to steeper SMB gradients and larger ELA variability.Since results suggest that areas with these characteristics are more sensitive to changes in SMB, we expect that modeling of outlet glaciers in Southern Greenland would require a resolution higher than 40 km.In the Southeast, the terrain is more mountainous, and outlet glaciers are narrow and steep.In this area, an even higher resolution in SMB may be required, as sensitivity results indicate that ice flux is very sensitive to the local driving stress.Indeed, we find that steeper areas (i.e., gates 4, 7, and 8) have larger mass flux uncertainties.Sensitivity results for gate 4 specifically offer insight into the behavior of ice flow through narrow channels.In Table 2, we find that the first-order radius of gate 4 is 16 km, a value much smaller than the other gates, suggesting that Southeast Greenland may indeed require finer resolution of SMB forcing, on the order of 15 km.These conclusions have strong implications for coupling between ice flow models and GCMs, as improved resolution and more accurate simulation of polar climate conditions are important to the advancement of ice flow model projections.
[44] Here we present results that heavily rely upon computationally demanding UQ methods.Considering the sizeable resources required by these analyses, SSA offers an advantage over more intensive approaches such as full-Stokes and Blatter-Pattyn models [Pattyn, 1996;Morlighem et al., 2010;Price et al., 2011;Larour et al., 2012c].However, it is important to acknowledge that in choosing the SSA, we make an assumption that vertical shear stress is negligible, which could potentially impact our results.Past ISSM ice stream studies have set the precedence for using SSA [Larour et al., 2012b[Larour et al., , 2012a]].Seroussi et al. [2011], for instance, evaluated Blatter-Pattyn and full-Stokes models of the northernmost outlet of the NEGIS, 79north Glacier.They determined that in most of the 79north domain, sliding is dominant, as the average difference between surface and depth-averaged velocity is on the order of 1%.Their results suggest that within the interior of the ice stream, where velocities are in excess of 60 m/yr, it is valid to assume that vertical shear stress is negligible.Over the entire NEGIS domain used in this study, we find that the majority of fast-flowing ice in the interior of the ice stream is indeed affected by sliding (Figures 2b and 2c), in agreement with similar modeling efforts in this region [Joughin et al., 2001].For our study domain, we carried out a comparison between the SSA model and a higher-order (Blatter-Pattyn) model.In areas where the basal drag coefficient has a value of 35 (m/s) -1/2 or less (Figure 2c), we determined that the average difference between surface and basal velocity is 1% of the depth-average velocity.This analysis suggests that SSA is a reasonable alternative to higher-order models in the simulation of high-velocity features such as the NEGIS.In agreement with Seroussi et al. [2011], we find that larger deviations occur in areas of slow moving ice along the margins.This affects the most narrow gates (3, 4, 7, and 8) along their profile margins, where there is an intense transition from high-velocity to low-velocity flow.We also find that SSA starts to break down in the upper part of the ice stream, where vertical shearing is more intense (surface to basal velocity differences reaching up to 6% of the depth-average velocity).Thus, we expect ice flux through upstream gates (gates 1, 2, and 5) to be less accurately simulated.These limitations should be taken into account when considering the results discussed above.

Conclusion
[45] In this study, we apply the newly developed ISSM framework to the Northeast Greenland Ice Stream and its Zachariae Isstrøm outlet.We extend the established uncertainty quantification capabilities, originally used to assess a static model, to a transient ice sheet model, by applying sampling and sensitivity analysis to a 22 year run forced by regional-climate-model-derived historical SMB.Our results demonstrate that errors in SMB affect mass flux in a way that is very similar to how errors in ice thicknesses propagate through a static inversion, suggesting that SMB effects on the driving stress dominate the model response.Extension to transient forcing does not expose instability issues in the model but instead further emphasizes how it serves as a filter that dampens input noise and smooths velocities, even after thousands of time steps.In addition, sensitivity analysis results illustrate that (1) local errors in SMB can affect ice flow up to 200 km away; (2) highly resolved and accurate surface forcing is critical in proximity to the ELA; and (3) in terms of coupling, ice flow simulation requires SMB forcing to be at a resolution of least 40 km.These conclusions have strong implications for coupling between ice flow models and GCMs, as relatively high-resolution and accurate representation of polar climate conditions are imperative for realistic ice flow simulation.

Figure 1 .
Figure1.Greenland Ice Sheet standard error (e) serving as input to the DAKOTA sampling analysis (m/yr water equivalent).In the ablation zone, e = 0.6 m/yr, while accumulation zone errors are equal to one standard deviation of the regression residual of accumulation reconstruction[Box et al., 2012;Box, 2013].Outlined in gray is the NEGIS regional domain (ISSM-NEGIS) of the 22 year run used in uncertainty quantification.This domain is extracted from the larger continental domain after completion of a 150 year relaxation.

Figure 2 .
Figure 2. (a) The model mesh, (b) modeled velocities at initialization of sampling and sampled mass flux gates, (c) basal drag coefficient, and (d) mean SMB over the sampling period (1989-2010).The state of the domain, extracted from the full Greenland domain after 150 years of relaxation, dictates the initial condition for the 22 year run assessed for uncertainty.

Figure 4 .
Figure 4. Sensitivity (Â, Gt/yr) of gates 1-4 (Figure 2b) relative to error in SMB.The domain is partitioned into 500 equal regions.SMB in each partition is perturbed by 0.1% and used as forcing for a distinct model run.Sensitivities (output:input) assess the impact of SMB errors on the mass flux through each gate.

Figure 5 .
Figure 5. Sensitivity (Â, Gt/yr) of gates 5-8 (Figure2b) relative to error in SMB.The domain is partitioned into 500 equal regions.SMB in each partition is perturbed by 0.1% and used as forcing for a distinct model run.Sensitivities (output:input) assess the impact of SMB errors on the mass flux through each gate.

Figure 6 .
Figure 6.Scaled sensitivities (SS) of gates 1-4 (Figure 2b), calculated by scaling the square of the sensitivity values plotted in Figure 4. Scaled sensitivities sum to unity over the domain.Results are displayed on a logarithmic scale in order to emphasize the extent of local SMB influence on each flux gate.

Figure 7 .
Figure 7. Scaled sensitivity (SS) of gates 5-8 (Figure 2b), calculated by scaling the square of the sensitivity values plotted in Figure 5. Scaled sensitivities sum to unity over the domain.Results are displayed on a logarithmic scale in order to emphasize the extent of local SMB influence on each flux gate.

Table 1 .
The Correlation Coefficient Between the Sampling Analysis Flux Responses at All Gates a

Table 2 .
The Maximum Radius of Influence (km) of SMB Forcing Errors on Ice Flux through Eight Gates, as Measured by Scaled Sensitivities (Figures6 and 7) a