Modelled glacier dynamics over the last quarter of a century at Jakobshavn Isbræ

. Observations over the past 2 decades show substantial ice loss associated with the speed-up of marine-terminating glaciers in Greenland. Here we use a regional three-dimensional outlet glacier model to simulate the behaviour of Jakobshavn Isbræ (JI) located in western Greenland. Our approach is to model and understand the recent behaviour of JI with a physical process-based model. Using atmospheric forcing and an ocean parametrization we tune our model to reproduce observed frontal changes of JI during 1990–2014. In our simulations, most of the JI retreat during 1990–2014 is driven by the ocean parametrization used and the glacier’s subsequent response, which is largely governed by bed geometry. In general, the study shows signif-icant progress in modelling the temporal variability of the ﬂow at JI. Our results suggest that the overall variability in modelled horizontal velocities is a response to variations in terminus position. The model simulates two major accelerations that are consistent with observations of changes in glacier terminus. The ﬁrst event occurred in 1998 and was triggered by a retreat of the front and moderate thinning of JI prior to 1998. The second event, which started in 2003 and peaked in the summer 2004, was triggered by the ﬁnal break-up of the ﬂoating tongue. This break-up reduced the buttressing at the JI terminus that resulted in further thinning. As the terminus retreated over a reverse bed slope into deeper water, sustained high velocities over the last decade have been observed at JI. Our model provides evidence that the 1998 and 2003 ﬂow accelerations are most likely initiated by the ocean parametrization used but JI’s subsequent dynamic response was governed by its own bed geometry. We are un-able to reproduce the observed 2010–2012 terminus retreat in our simulations. We attribute this limitation to either inaccuracies in basal topography or to misrepresentations of the climatic forcings that were applied. Nevertheless, the model is able to simulate the previously observed increase in mass loss through 2014.


Introduction
The rate of net ice mass loss from Greenland's marineterminating glaciers has more than doubled over the past 2 decades (Rignot et al., 2008;Moon et al., 2012;Shepherd et al., 2012;Enderlin et al., 2014). Jakobshavn Isbrae, located mid-way up the western side of Greenland, is one of the largest outlet glaciers in terms of drainage area as it drains ∼ 6 % of the Greenland Ice Sheet (GrIS) (Krabill et al., 2000). Due to its consistently high ice flow rate and seasonally varying flow speed and front position, the glacier has received much attention over the last 2 decades (Thomas et al., 2003;Luckman and Murray, 2005;Holland et al., 2008;Amundson et al., 2010;Khan et al., 2010; Joughin et al., 2012;Gladish et al., 2015a, b;de Juan et al., 2010). Measurements from synthetic aperture radar suggest that the ice flow speed of JI doubled between 1992 and 2003 (Joughin et al., 2004). More recent measurements show a steady increase in the flow rate over the glacier's fastermoving region of ∼ 5 % per year (Joughin et al., 2008). The speed-up coincides with thinning of up to 15 m a −1 between 2003 and 2012 near the glacier front (Krabill et al., 2004;Nielsen et al., 2013) as observed from airborne laser altimeter surveys. The steady increase in the flow rate and glacier thinning suggests a continuous dynamic drawdown of mass and highlights JI's importance for the GrIS mass balance.
Over the past decade, we have seen significant improvements in the numerical modelling of glaciers and ice sheets (e.g. Price et al., 2011;Vieli and Nick, 2011;Larour et al., 2012;Pattyn et al., 2012;Seroussi et al., 2012;Aschwanden et al., 2013Aschwanden et al., , 2016Nick et al., 2013;Mengel and Levermann, 2014) and several processes have been identified as controlling the observed speed-up of JI (Nick et al., 2009;Van der Veen et al., 2011;. One process is a reduction in resistance (buttressing) at the marine front through thinning and/or retreat of the glacier termini. However, the details of the processes triggering and controlling thinning and retreat remain elusive. Accurately modelling complex interactions between thinning, retreat, and acceleration of flow speed as observed at JI is challenging. Our knowledge of the mechanisms triggering these events is usually constrained to the period covered by observations. The initial speed-up of JI occurred at a time when the satellite and airborne observations were infrequent and therefore insufficient to monitor the annual to seasonal evolution of glacier geometry and speed.
Here, we use a high-resolution, three-dimensional, timedependent regional outlet glacier model that has been developed as part of the Parallel Ice Sheet Model (PISM; see Sect. 2.1) (Khroulev and the PISM Authors, 2014) to investigate the dynamic evolution of JI between 1990 and 2014. While previous three-dimensional modelling studies have mostly concentrated on modelling individual processes using stress perturbations (e.g. Van der Veen et al., 2011;Joughin et al., 2012), the present study aims to model the recent behaviour of JI with a process-based model. Our modelling approach is based on a regional equilibrium simulation and a time integration over the period 1990 to 2014, in which the grounding lines and the calving fronts are free to evolve under the applied ocean parametrization and monthly atmospheric forcing.

Ice sheet model
The ice sheet model used in this study is the PISM (stable version 0.6). PISM is an open-source, parallel, three-dimensional, thermodynamically coupled, and timedependent ice sheet model (Bueler and Brown, 2009;Khroulev and the PISM Authors, 2014). The model uses the superposition of the non-sliding shallow ice approximation (SIA; Hutter, 1983) for simulating slowly moving grounded ice in the interior part of the ice sheet and the shallow shelf approximation (SSA; Weis et al., 1999) for simulating fastflowing outlet glaciers and ice shelf systems. We solve the SIA with a non-sliding base and use the SSA as a basal sliding velocity for the ice grounded regions . This superposition of SIA and SSA (the "SIA + SSA" hybrid model) sustains a smooth transition between nonsliding, bedrock-frozen ice and sliding, fast-flowing ice and has been shown to reasonably simulate the flow of both grounded and floating ice . To determine driving stresses for the SIA and SSA stress balances, PISM computes surface gradients according to Mahaffy (1976). For conservation of energy, we use an enthalpy scheme (Aschwanden et al., 2012) that accounts for changes in temperature in cold ice (i.e. ice below the pressure melting point) and for changes in water content in temperate ice (i.e. ice at the pressure melting point).
In PISM, the basal shear stress is related to the sliding velocity through a nearly plastic power law (Schoof and Hindmarsh, 2010). The Mohr-Coulomb criterion (Cuffey and Paterson, 2010) is used to connect a saturated and pressurized subglacial till with a modelled distribution of yield stress. The yield stress depends on the effective pressure and on a spatially varying till friction angle derived heuristically as a piecewise-linear function of the bed elevation Winkelmann et al., 2011;Aschwanden et al., 2013). The effective pressure on the till is determined by the ice overburden pressure and the effective thickness of water in the till (Tulaczyk et al., 2000a, b). In this subglacial hydrology model the water is not conserved and it is only stored locally in the till up to a maximum thickness of 2 m. The ice flow therefore develops in PISM as a consequence of plastic till failure, i.e. where the basal shear stress exceeds the yield stress, and is influenced by the thermal regime and the volume of water at the ice sheet bed.
The underlying equations are further illustrated in the Supplement.

Input data
We use the bed topography from Bamber et al. (2013). This 1 km bed elevation data set for all of Greenland was derived from a combination of multiple airborne ice thickness surveys and satellite-derived elevations during 1970-2012. The data set has an increased resolution along the ice sheet margin. In the region close to the outlet of JI, data from an 125 m CReSIS digital elevation model (that includes all the data collected in the region by CReSIS between 1997 and2007) have been used to improve the accuracy of the data set. Errors in bed elevation range from 10 to 300 m, depending on The Cryosphere, 10, 597-611, 2016 www.the-cryosphere.net/10/597/2016/ the distance from an observation and the variability of the local topography . The terminus position and surface elevation in the Jakobshavn region are based on 1985 aerial photographs (Csatho et al., 2008). Ice thickness in the JI basin is computed as the difference between surface and bedrock elevation. The model of the geothermal flux is adopted from Shapiro and Ritzwoller (2004). We use monthly input fields of near-surface air temperature and surface mass balance (SMB) from the regional climate model RACMO2.3 (Noël et al., 2015;Figs. S2 and S3 in the Supplement), which here represent the only seasonal input used in the model. The version used in this study is produced at a spatial resolution of ∼ 11 km and covers the period from 1958 to 2014. Additional grid refinements are performed using bilinear interpolation for climatic data sets and a secondorder conservative remapping scheme (Jones, 1999) for bed topography data.

Initialization procedure, boundary conditions, calving, and grounding line parametrization
In our model, the three-dimensional ice enthalpy field, basal melt for grounded ice, modelled amount of till-pore water, and lithospheric temperature are obtained from an ice-sheetwide palaeoclimatic spin-up. The palaeoclimatic spin-up follows the initialization procedure described by Bindschadler et al. (2013) and Aschwanden et al. (2013). We start the spinup on a 10 km grid, and then we further refine it to 5 km at −5 ka. It is important to note that during the palaeoclimatic initialization the terminus is held fixed to the observed 1990 position in the JI region and to the position from  elsewhere. In the regional outlet glacier model of PISM, the boundary conditions are handled in a 10 km strip positioned outside of the JI's drainage basin and around the edge of the computational domain (Fig. 1b). In this strip, the input values of the basal melt, the amount of till-pore water, ice enthalpy, and lithospheric temperature  are held fixed and applied as Dirichlet boundary conditions in the conservation of energy model (Khroulev and the PISM Authors, 2014). The boundary conditions for the enthalpy at the ice-bedrock interface follow Aschwanden et al. (2012). We start our regional JI runs with an equilibrium simulation on a horizontal grid with 5 km spacing. The enthalpy formulation models the mass and energy balance for the threedimensional ice fluid field based on 200 regularly spaced ice layers within a domain extending 4000 m above the bed elevation. The temperature of the bedrock thermal layer is computed up to a depth of 1000 m with 50 regularly spaced layers. The first step is to obtain a 5 km regional equilibrium model for JI using constant mean climate (i.e. repeating the 1960-1990 mean air temperature and SMB; see Sect. 2.1.1). We consider that equilibrium has been established when the ice volume in the regional domain changes by less than 1 % in the final 100 model years. Grid refinements are made from 5 km (125 × 86) to 2 km (310 × 213) after 3000 years. The 2 km simulation reaches equilibrium after 200 years with an ice volume of 0.25 × 10 6 km 3 (or a 3.6 % increase relative to the input data set from Bamber et al., 2013). Further, using our equilibrium simulations with a 2 km horizontal grid and 400 regularly spaced ice layers within a domain extending 4000 m above the bed elevation, we simulate forward in time (hindcast) from 1990 to 2014 by imposing monthly fields of SMB and 2 m air temperatures through a one-way forcing scheme. For simulations performed on a 1 km horizontal grid, the exact same procedure is used with the additional constraint that in the regional equilibrium run a further grid refinement from 2 to 1 km is made after 200 years. The length of the 1 km regional equilibrium simulation is 100 years. In our regional model, all boundaries (calving fronts, grounding lines, upper, and lower surfaces) are free to evolve in time both during the regional equilibrium and the forward simulations. Along the ice shelf calving front, we superimpose a physically based calving (eigen-calving) parametrization Levermann et al., 2012) and www.the-cryosphere.net/10/597/2016/ The Cryosphere, 10, 597-611, 2016 a basic calving mechanism ) that removes any floating ice at the calving front thinner than a given threshold at a maximum rate of one grid cell per time step. The average calving rate (c) is calculated as the product of the principal components of the horizontal strain rates (ε ± ), derived from the SSA velocities, and a proportionality constant parameter (k) that captures the material properties relevant for calving: The strain rate pattern is strongly influenced by the geometry and the boundary conditions at the ice shelf front (Levermann et al., 2012). The proportionality constant, k, is chosen such that the ice front variability is small (Leverman et al., 2012). This physically based calving law appears to yield realistic calving front positions for various types of ice shelves having been successfully used for modelling calving front positions in entire Antarctica simulations  and regional east Antarctica simulations (Mengel and Levermann, 2014). In contrast to Antarctica, known for its large shelves and shallow fjords, the GrIS is characterized by narrow and deep fjords, and JI is no exception. The strain rate pattern in the eigen-calving parametrization performs well only if fractures in glacier ice can grow, and calving occurs only if these rifts intersect (i.e. possible only for relatively thin and unconfined ice shelves). In the case of JI, whose terminus is confined in a narrow fjord, the strain rate pattern that defines the eigen-calving parametrization is not the governing process. In our model, the eigen-calving law has priority over the basic calving mechanism. That is to say that the second calving law used (the basic calving mechanism) removes any ice at the calving front not calved by the eigen-calving parametrization thinner than 500 m in the equilibrium simulations and 375 m in the forward runs. Therefore, the creation of the conditions under which calving can occur (e.g. a floating ice shelf) with the subsequent calving mechanism relies solely on the parametrization for ice shelf melting (Sect. 2.1.3).
A partially filled grid cell formulation , which allows for sub-grid-scale retreat and advance of the ice shelf front, is used to connect the calving rate computed by the calving parametrizations with the mass transport scheme at the ice shelf terminus. This sub-grid-scale retreat and advance of the shelf allows for realistic spreading rates that are important for the eigen-calving parametrization. The sub-grid interpolation is performed only when a floating terminus exists. In both situations (i.e. floating ice or grounded terminus), the stress boundary conditions are applied at the calving front and in the discretization of the SSA equations . The retreat and advance of the front through calving is restricted to at most one grid cell length per adaptive time step.
The parametrization of the grounding line position is based on a linear interpolation scheme (the "LI" parametrization; Gladstone et al., 2010) extended to two horizontal dimensions (xy) and is not subject to any boundary conditions. This sub-grid treatment of the grounding line interpolates the basal shear stress in xy based on the spatial gradient between cells below and above the grounding line and allows for a smooth transition of the basal friction from grounded to floating ice . At each time step the grounding line position is determined by a mask that distinguishes between grounded and floating ice using a flotation criterion based on the modelled ice thickness where b represents the bedrock elevation, ρ i is the density of the ice, ρ o is the density of the ocean water, and H represents the ice thickness. Therefore, the grounding line migration is influenced by the ice thickness evolution, which further depends on the velocities computed from the stress balance. The superposition of SIA and SSA, which implies that the SSA velocities are computed simultaneously for the shelf and for the sheet, ensures that the stress transmission across the grounding line is continuous and that buttressing effects are included. In the three-dimensional Marine Ice Sheet Model Intercomparison Project (Mismip3d), PISM was used to model reversible grounding line dynamics and produced results consistent with full Stokes models ; see parameters therein).
We have not performed the Mismip3d experiments for our particular parameter settings and, therefore, the accuracy of the modelled grounding line migration is solely based on the results presented in .

Parametrization for ice shelf melting
We use a simple parametrization for ice shelf melting where the melting effect of the ocean is based on both sub-shelf ocean temperature and salinity . To accommodate this parametrization, several changes have been made to PISM at the sub-shelf boundary . First, the ice temperature at the base of the shelf (the pressure-melting temperature) necessary for the enthalpy solver (Aschwanden et al., 2012) is calculated from the Clausius-Clapeyron gradient and the elevation at the base of the shelf. The ice temperature is then applied as a Dirichlet boundary condition in the conservation of energy equation. Secondly, basal melting and refreezing is incorporated through a sub-shelf mass flux used as a sink/source term in the mass-continuity equation. This mass flux from shelf to ocean (Beckmann and Goosse, 2003) is computed as a heat flux between the ocean and ice and represents the melting effect of the ocean due to both temperature and salinity .
In our simulations we use a constant ocean water temperature (T o ) of −1.7 • C, which here represents the mean sur-The Cryosphere, 10, 597-611, 2016 www.the-cryosphere.net/10/597/2016/ face ocean temperature in the grid cells adjacent to the JI terminus. In the heat flux parametrization, the ocean temperature at the ice shelf base is computed as the difference between the input ocean temperature and a virtual temperature that represents the freezing point temperature of ocean water below the ice shelf (Fig. S4). The freezing point temperature is calculated based on the elevation at the base of the shelf and the ocean water salinity. As a consequence of these constraints, as the glacier retreats and/or advances both the pressure-melting temperature and the heat flux between the ocean and ice evolve alongside the modelled glacier ice shelf geometry. The ocean water salinity (S o = 35 psu) is kept constant in time and space as the model does not capture the salinity gradient from the base of the ice shelf through layers of low and high salinity. A previous study conducted by Mengel and Levermann (2014) using the same model established that the sensitivity of the melt rate to salinity is negligible. Following for this melting parametrization, the highest melt rates are modelled in the proximity of the glacier grounding lines and decrease with elevation such that the lowest melt rates are closer to the central to frontal area of the modelled ice shelf. At the grounding line, PISM computes an extra flotation mask that accounts for the fraction of the cell that is grounded by assigning 0 to cells with fully grounded ice, 1 to cells with ice-free or fully floating ice, and values between 0 and 1 to partially grounded grid cells. The basal melt rate in the cells containing the grounding line is then adjusted based on this flotation mask as follows (Khroulev and the PISM Authors, 2014): where M b refers to the basal melt rate and λ is the value of the flotation mask. At the vertical ice front, we do not apply any melt.

Results and discussion
This section is organized in two main subsections. Section 3.1 introduces the results obtained relative to observations, and Sect. 3.2 focuses mainly on the limitations of the model that need to be considered before a final conclusion can be drawn. A short introduction to the different simulations and preparatory experiments performed is given below. A total number of 50 simulations with different sets of parameters (excluding preparatory and additional experiments on the 1 km grid) are performed on a 2 km grid. We alter six parameters that control the ice dynamics (e.g. the flow enhancement factor, the exponent of the pseudo-plastic basal resistance model, the till effective fraction overburden), the ice shelf melt, the ocean temperature, and the calving (i.e. the ice thickness threshold in the basic calving mechanism). These parameters are modified only during the regional JI runs such that the model reproduces the frontal positions and the ice mass change observations at JI during the period From these results, we present the parametrization that best captures (i.e. we estimate the residual between modelled and observed ice mass change and select the smallest residual signal) the full observed evolution of JI during the period 1990-2014 (Figs. 2, 3, 4, and5). The values of the ice sheet model parameters used and the ice sheet model sensitivity to parameters controlling ice dynamics, basal processes, ice shelf melt, and ocean temperature are further illustrated in the Supplement.

Annual-scale variations in velocities, terminus, and grounding line positions
We investigate the processes driving the dynamic evolution of JI and its variation in velocity between 1990 and 2014 with a focus on the initial speed-up of JI (1990) figure) for the period 1990-2014 at locations (S1 to S7) shown in Fig. 1c. The same colour scheme is used for the modelled and the observed data. The observed velocities prior to 2009 are mean winter velocities and are largely consistent with our modelled winter estimates for the same period. The observed thickness has been adjusted to match the model thickness at the first available observation (i.e. by summing the modelled ice thickness corresponding to the first available observation with the observed thickness changes).

1990-1997
The first speed-up produced by the simulation is caused by a retreat of the front position by approximately 2 to 4 km between 1990 and 1991. There is no observational evidence to confirm that this retreat actually occurred. The simulated retreat is probably a modelling artefact as the geometry obtained during the regional equilibrium simulation is forced with monthly atmospheric forcing and new oceanic conditions. This simulated acceleration (Fig. 3) is caused in our model by a reduction in buttressing due to a reduction in lateral resistance (Van der Veen et al., 2011), which is generated by the gradual retreat of the front and triggers a dynamic response in the upstream region of JI. Starting in 1992, the modelled and observed terminus positions agree (not shown in Fig. 2). Apart from the acceleration in 1991-1992, no significant seasonal fluctuations in flow rate are found in our simulations for this period, a result that is consistent with observations (Echelmeyer et al., 1994). From 1993 a stronger sub-annual velocity signal begins to emerge in our simulation that continues and intensifies in magnitude during 1994 and 1995. Modelled mean-annual ve-locities for 1992 and 1995 are consistent with observed velocities for the same period (Joughin et al., 2008;Vieli and Nick, 2011). In 1996 and 1997, the frontal extent and the grounding line position remain relatively stable (Figs. 2, 6, and 7), and no significant seasonal fluctuation in ice flow rate is observed in the simulation. These model results agree well with observations, which indicate that the glacier speed was relatively constant during this period (Luckman and Murray, 2005).

1998-2002
According to observations (Joughin et al., 2004;Luckman and Murray, 2005;Motyka et al., 2011;Bevan et al., 2012), the initial acceleration of JI occurred in May-August 1998, which coincides with our modelled results. In our simulation, the 1998 acceleration is generated by a retreat of the ice tongue's terminus in 1997-1998, which may be responsible for reducing buttressing (Fig. 7 and Movie 1 in the Supplement). Thinning, both near the terminus and inland (up to 10 km away from the 1990 front position), starts in our model in the summer of 1995 and continues to accelerate after 1998 (Figs. 3, 6, and 7). The modelled behaviour agrees The Cryosphere, 10, 597-611, 2016 www.the-cryosphere.net/10/597/2016/ well with the observed behaviour (Krabill et al., 2004). Although thinning appears to have increased in our model during 3 continuous years, it produced only minor additional speed-up during the period prior to 1998 (Figs. 2, 6, and 7). In our simulation, JI's speed increased in the summer of 1998 by ∼ 80 % relative to the summer of 1992 (Fig. 3), at which time the grounding line position starts to retreat thereafter (Figs. 2, 6, and 7). Observations (Luckman and Murray, 2005) do not show this level of speed-up, and there are no observations of the grounding line position at this time with which to assess our model performance. Overall, modelling results suggest an advance of the terminus between 1999 and 2000 and a retreat of the southern tributary between 2000 and 2002 by ∼ 4 km, which correlates with existing observations (Thomas, 2004). In our simulation, this retreat of the terminus triggers a decrease of resistive stresses at the terminus (Figs. 7 and S8). Concurrent with the 1998-2002 terminus retreat, the grounding line retreats in our model by ∼ 6 km (Figs. 2, 6, and 7).

2003-2014
In the late summer of 2003, the simulated flow velocity increases (Fig. 3). This acceleration of JI is driven in  (Figs. 3 and 6) both near the front and further inland in response to a change in the near-terminus stress field (Fig. 7). During the final break-up of the ice tongue, the simulation produces speeds high as 20 km a −1 (∼ 120 % increase relative to 1998). The modelled velocities decreased to 16 km a −1 (∼ 80 % increase relative to 1998) in the subsequent months and remained substantially higher than the sparse observations from that time (e.g. Joughin et al., 2012). The high velocities modelled at JI after the loss of its floating tongue are further sustained in our simulation by the thinning that occurred after 2003 (Fig. 3), which continues to steepen the slopes near the terminus (Fig. 6). This simulated thinning is accompanied by a seasonal driven (sub-annual scale) retreat and advance of the front and is combined in the following years with a reduction in surface mass balance due to in- is characterized in our simulation by relatively uniform velocity peaks with strong sub-annual variations (Fig. 3). During this period, only a small floating ice tongue is modelled and the terminus remained relatively stable, with no episodes of significant retreat.
In agreement with previous studies (e.g. Joughin et al., 2012), our results suggest that the overall variability in the modelled horizontal velocities is a response to variations in terminus position (Fig. 7). In our simulation, the retreat of the front reduced the buttressing at the terminus and generated a dynamic response in the upstream region of JI which finally led to flow acceleration. In contrast, when the front advanced the modelled flow slowed as the resistive stresses at the terminus were reinforced. This buttressing effect tends to govern JI's behaviour in our model. Regarding the overall terminus retreat, our simulations suggest that it is mostly driven by the sub-shelf melting parametrization applied (Figs. S5 and S14). Although the heat flux supplied to the shelf evolves in time based on the modelled terminus geometry, the input ocean temperature is kept constant throughout the simulations. This constant ocean forcing at the terminus leads, in our simulation, to gradual thinning of JI and favours its retreat without any shift (e.g. increase) in ocean temperature. In terms of seasonality, the only seasonal input into the model is introduced by the monthly atmospheric forcing that is applied (Sect. 2.1.1). In our model, the atmospheric forcing that is applied (Figs. S2 and S3) can influence JI's dynamics through changes in SMB (i.e. accumulation and ablation), which affects both the SIA and the SSA (Sect. 2.1). However, the modelled sub-annual variability in terms of terminus retreat and velocities does not always follow the seasonal signal (Fig. 3). We investigate this higher than seasonal variability in Sect  Figure 4 shows observed and modelled mass change for the period 1997 to 2014. We estimate the observed rate of ice volume changes from airborne and satellite altimetry over the same period and convert these to rates of mass change (Supplement, Sect. 2). Overall we find good agreement between modelled and observed mass change (Fig. 4), and our results are in agreement with other similar studies Nick et al., 2013). Dynamically driven discharge is known to control Jakobshavn's mass loss between 2000 and 2010 (Nick et al., 2013). The modelled cumulative mass loss is 269 Gt, of which 93 % (∼ 251 Gt) is dynamic in origin while the remaining 7 % (∼ 18 Gt) is attributed to a decrease in SMB (Fig. 4). Further, the present-day unloading of ice causes the Earth to respond elastically. Thus, we can use modelled mass changes to predict elastic uplift. We compare modelled changes of the Earth's elastic response to changes in ice mass to uplift observed at four GPS sites (Fig. 5). Both model predictions and observations consistently suggest large uplift rates near the JI front (20 mm a 1 for station KAGA) and somewhat minor uplift rates (∼ 5 mm a −1 ) at distances of > 100 km from the ice margin.

Ice mass change
Although the terminus has ceased to retreat in our simulations after 2009 (Figs. 6 and 7), the modelled mass loss, and more importantly the dynamic mass loss, continues to accelerate (Fig. 4). Our results show (Fig. 7) that during this period the mass change is mostly driven by the sub-annual terminus retreat and advance, which continues to generate dynamic changes at JI through seasonal (sub-annual scale) reductions in resistive stresses.

Feedback mechanisms, forcings, and limitations
Representing the processes that act at the marine boundary (i.e. calving and ocean melt) is important for understanding and modelling the retreat/advance of marine-terminating glaciers like JI. Determining terminus positions by using the superposition of a physically based calving (eigen-calving) parametrization Levermann et al., 2012) and a basic calving mechanism  is motivated by the model's ability to maintain realistic calving front positions (Levermann et al., 2012). The eigencalving parametrization cannot resolve individual calving events, and, thus, the introduction of the basic calving mechanism was necessary in order to accurately match observed front positions. Preparatory experiments have shown that calving is mostly driven in our model by the basic calving mechanism used (∼ 96 % of the overall mass loss) and that the eigen-calving parametrization is more important in modelling sub-annual to seasonal fluctuations of the terminus. Our simulations suggest that the superposition of these two calving mechanisms performs well for relatively narrow and deep fjords as those characterized by JI (Fig. 2). The benefit of using such a combination of calving laws is that it can evolve the terminus position with time and thus calving feedbacks are not ignored. As the terminus retreats, the feedback between calving and retreat generates dynamic changes due to a reduction in lateral shear and resistive stresses (Fig. 7). In a simulation in which the terminus position is kept fixed to the 1990s position, the velocity peaks are uniform (i.e. no acceleration is modelled except for some small seasonal related fluctuations generated by the atmospheric forcing applied), and the mass loss remains relatively small (∼ 70 Gt). , 10, 597-611, 2016 www.the-cryosphere.net/10/597/2016/ Consistent with Vieli et al. (2011), we find that the feedback between calving and retreat is highly important in modelling JI's dynamics. As introduced in Sect. 2, our approach here is to adjust the terminus in the JI region to simulate the 1990s observed front position and surface elevation based on 1985 aerial photographs (Csatho et al., 2008). The glacier terminus in 1990s was floating (Csatho et al., 2008;Motyka et al., 2011). Motyka et al. (2011) calculated the 1985 hydrostatic equilibrium thickness of the south branch floating tongue from smoothed surface digital elevation models and obtained a height of 600 m near the calving front and 940 m near the grounding zone. In this paper, however, we compute the thickness as the difference between the surface elevation and the bed topography and allow the glacier to evolve its own terminus geometry during the equilibrium simulation. Preparatory experiments have shown that in our model (disregarding its initial geometry floating/grounded terminus) JI attains equilibrium with a grounding line position that stabilizes close to the 1990s observed terminus position. According to observations, JI is characterized in 1990 by a large floating tongue (> 10 km; e.g. Motyka et al., 2011) that we are not able to simulate during the equilibrium runs. In our model (Figs. 6  and 7), the glacier starts to develop a large floating tongue (∼ 10 km) in 1999. Starting in 2000, the floating tongue is comparable in length and thickness with observations and the model is able to simulate, with a high degree of accuracy, its break-up that occurred in late summer 2003 and the subsequent glacier acceleration. Observations of terminus positions (Sohn et al., 1998;Csatho et al., 2008) suggest that over more than 40 years, between 1946 and 1992, JI's terminus stabilized in the proximity of the 1990's observed terminus position. Furthermore, during 1959 and 1985 the southern tributary was in balance (Csatho et al., 2008). This suggests that, during the regional equilibrium and at the beginning of the forward simulations, we are forcing our model with climatic conditions that favoured the glacier to remain in balance. This may explain our unsuccessful attempts to simulate prior to 1998 a floating tongue comparable in length and thickness with observations and suggests that future studies should consider modelling JI before the glacier begins to float in the late 1940s to simulate the large floating tongue that characterized JI during this period (Csatho et al., 2008).

The Cryosphere
The geometry of the terminus plays an important role in parametrizing ice shelf melting, and therefore our pre-1999 geometry will influence the magnitude of the basal melt rates (Sect. 2.1.3). The difference in geometry results in modelled mean basal melt rates that are larger for the period 1999-2003 (Table S3), when JI begins to develop a large floating tongue and when the calving front was already largely floating. The modelled mean melt rates for the period 1999-2003 are large and likely overestimated. Relative to other studies, e.g. Motyka et al. (2011), our yearly mean melt rate for 1998 is ∼ 2 times larger (Table S3). While we choose here to compare the two melt rates in order to offer a scale perspective, we acknowledge the difference in geometry between the two studies.
Starting in 2010, the retreat of the terminus modelled in our simulations did not correlate well with observations (Fig. 2). The observed terminus and the grounding line retreats do not cease after 2010. Further, observed front positions (Joughin et al., 2014)  Modelled horizontal velocities and ice thickness changes at the point location S1 shown in Fig. 1c. (c) Modelled two-dimensional deviatoric stresses (in the X direction, the Y direction, and the shear stress) at the point location S1 shown in Fig. 1c.
JI was already retreating over the sill and on the over deepening indicated by the red star in Fig. 6. The observed retreat is not reproduced in our simulations suggesting that additional feedbacks and/or forcings most likely affect the glacier. Alternatively, the mismatch between observations and simulation results may represent an incomplete modelling of the physics, inaccuracies in atmospheric/oceanic conditions, or other various limitations (e.g. bed topography model constraints and grid resolution constraints). The particular influence of these potential limitations on our model is detailed below. The basal topography of JI's channels represents a large source of uncertainty. JI is a marine-terminating glacier whose bedrock topography is characterized by a long and narrow channel with deep troughs that contribute to its retreat and acceleration; e.g. once the grounding line starts to retreat on a down-sloping bed the flow increases, leading to further retreat and acceleration (Vieli et al., 2011). The timing and the magnitude of these retreats depend on bed topography and the glacier width changes (Jamieson et al., 2012;Enderlin et al., 2013). Accurate modelling of the grounding line behaviour is, therefore, crucial for JI's dynamics as its retreat removes areas of flow resistance at the base and may trigger unstable retreat if the glacier is retreating into deeper waters. In our simulation, the grounding line position stabilizes downstream of the sill after 2005 (Figs. 2 and 6), which is in accordance with previous modelling studies (Vieli et al., 2001;Vieli and Nick, 2011). Vieli and Nick (2011) found that, by artificially lowering the same bed sill by 100 m, the grounding line eventually retreats and triggers a catastrophic retreat of 80 km in just over 20 years. In an equivalent experiment with Vieli and Nick (2011) but performed with our model, lowering the bed sill by 100 m did not result in a retreat of the grounding line over the sill. Regarding the grid resolution, simulations performed on a 1 km grid did not improve our simulations of ice thickness (Fig. S10) or surface speed (i.e. trend, overall magnitude, and shape of the flow; Fig. S11).
From a climatic perspective, the summer of 2012 was characterized by exceptional surface melt covering 98 % of the entire ice sheet surface and including the high-elevation summit region (Nghiem et al., 2012;Hanna et al., 2014). Overall, the 2012 melt season was 2 months longer than the 1979-2011 mean and the longest recorded in the satellite era (Tedesco et al., 2013). Furthermore, the summer of 2012 was preceded by a series of warm summers (2007, 2008, 2010, and 2011) (Hanna et al., 2014). Surface melt above average was already recorded in May-June 2012 (see Fig. 3

from
The Cryosphere, 10, 597-611, 2016 www.the-cryosphere.net/10/597/2016/ NSIDC, 2015) when most of the 2011-2012 winter accumulation melted and over 30 % of the ice sheet surface experienced surface melt. An intense and long melt year leads to extensive thinning of the ice and has the potential to enhance hydrofracturing of the calving front due to melt water draining into surface crevasses (MacAyeal et al., 2003;Joughin et al., 2013;Pollard et al., 2015), resulting in greater and/or faster seasonal retreat and an increase in submarine melt at the terminus and the sub-shelf cavity (Schoof, 2007;Stanley et al., 2011;Kimura et al., 2014;Slater et al., 2015). The seasonal retreat of JI's terminus started relatively early in 2012, with a large calving event having already occurred in June. While it seems difficult to attribute this particular calving event solely to processes related to the 2012 melt season, it does seem probable that the series of warm summers (2007)(2008)(2009)(2010)(2011) together with the 2012 exceptional melt season could have enhanced hydrofracturing of the calving front. In turn, this could have induced a retreat of the terminus that cannot be captured by our model (i.e. in its present configuration the model cannot account directly for the influence of meltwater runoff and its role in the subglacial system during surface melt events). However, changes in ice thickness affect both the SIA and the SSA (Sect. 2.1). While the effect on the SIA is very weak as the driving stresses are not affected by a few metres of difference in thickness induced by SMB variability, in the SSA, the coupling is achieved via the effective pressure term in the definition of the yield stress (see Supplement, Sect. 1.2, for detailed equations). The effective pressure is determined by the ice overburden pressure (i.e. ice thickness) and the effective thickness of water in the till, where the latter is computed by time integrating the basal melt rate. Compared with SIA, this effect is stronger and may explain why in our model some seasonal velocity peaks could potentially be influenced by the atmospheric forcing applied (Figs. S9 and S14).
We study the sensitivity of the model to atmospheric forcing by performing a simulation where we keep the atmospheric forcing constant (mean 1960-1990 temperature andSMB). By comparing this simulation with a simulation that includes full atmospheric variability (monthly temperature and SMB) we find that to only a relatively small degree some of the variability appears to be influenced by the atmospheric forcing applied (Figs. S2 and S14), which also represents the only seasonal input into the model. Some of the greater than seasonal frequency could be an issue with resolution in the model. We examined this sensitivity by performing additional runs at a higher spatial resolution. Simulations on a 1 km grid did show some improvement with respect to surface speed sub-annual variability (Fig. S12), suggesting that in our model the stress redistribution might be sensitive to the resolution of the calving event. However, given the short period spanned by the simulations, the stress redistribution does not change the overall modelled results, as seen in Figs. S10 and S11. Although we acknowledge that some of the variability is due to the grid resolution, part of it may also be related to unmodeled physical processes acting at the terminus. We suggest that additional contributions to the seasonality, e.g. from ice mélange or seasonal ocean temperature variability, which are not included in our model, could potentially influence the advance and retreat of the front at seasonal scales (Fig. S14). For example, the ice mélange can prevent the ice at the calving front from breaking off and could therefore reduce the calving rates. Consequently, the introduction of an ice mélange parametrization will probably help to minimize some of the sub-annual signal modelled in our simulations. Similarly, seasonal ocean temperature variability can influence ice mélange formation and/or clearance and the melt rates at the glacier front and can accentuate seasonal glacier terminus and grounding line retreat and/or advance. However, at this point we find it difficult to determine the relative importance of each process.
Finally, regarding the ocean conditions, warm water temperatures in the fjord were recorded in 2012. Besides a cold anomaly in 2010, which was sustained until early 2011, the period 2008-2013 is characterized by high fjord water temperatures -equal to or warmer than those recorded in 1998-1999 (Gladish et al., 2015a, b). In our model, the ice melt rates are determined from the given conditions in temperature (−1.7 • C and salinity (35 psu) of the fjord waters) and the given geometry (Sect. 2.1.3). The fact that we are able to model JI's retreat with a constant ocean temperature suggests that the retreat and acceleration observed at JI are not likely to be controlled by the year-to-year variability in ocean temperatures. This conclusion agrees with the observational study of Gladish et al. (2015a, b), who analysed ocean temperature variability in the Ilulissat fjord with JI variability and found that after 1999 there was no clear correlation. Our results do not, however, imply that the ocean influence in JI's retreat is negligible (Fig. S5), but rather that the glacier most likely responds to changes in ocean temperature that are sustained for longer time periods, e.g. decadal timescales. Two additional experiments, in which the input ocean temperature (T o ) was increased to −1 • C indicate that higher melt rates beneath the grounding line could potentially explain the retreat observed after 2010. In our first experiment, the input T o was increased from −1.7 to −1 • C starting in 1997 (∼ 0.7 • C relative to 1990). This temperature increase is consistent with observed ocean temperatures at the mouth of the Ilulissat fjord (Gladish et al., 2015a, b) and generated in our simulation, for the period 1997-2014, an accelerated retreat of the front that does not correlate with observations (Fig. S7). Similarly, mass loss estimates from the simulations are significantly larger (by ∼ 50 %; Fig. S6) than those calculated from airborne and satellite altimetry observations (Sect. 3.1.2). Overall, the experiment shows that an increase in ocean temperature that starts in 1997 and is sustained until 2014 generates modelled estimates for the period 1998-2014 that do not agree with observations. In the second experiment, T o was increased to −1 • C starting in 2010 (∼ + 0. our model predicted a faster retreat of the front that correlates well with observations ( Fig. S7) and an increase of mass loss by ∼ 7 Gt (Fig. S6). This experiment shows that an increase in ocean temperature beginning in 2010 could potentially explain the retreat observed thereafter.

Conclusions
In this study, a three-dimensional, time-dependent regional outlet glacier model is used to investigate the processes driving the dynamic evolution of JI and its seasonal variation in ice velocity between 1990 and 2014. Here, we attempted to simulate the recent behaviour of JI with a process-based model. The model parameters were calibrated such that the model reproduced observed front positions (Fig. 2) and ice mass change observations (Fig. 4) at JI over the periods 1990-2014 and 1997-2014 respectively. We obtain a good agreement of our model output with time series of measured horizontal velocities, observed thickness changes, and GPSderived elastic uplift of the crust (Figs. 3 and 5). Overall, the study shows progress in modelling the temporal variability of the flow at JI. Our results suggest that most of the JI retreat during 1990-2014 is driven by the ocean parametrization and the glacier's subsequent response, which is largely governed by its own bed geometry (Figs. 6,7,and S5). In agreement with previous studies (e.g. Joughin et al., 2012), our simulations suggest that the overall variability in the modelled horizontal velocities is a response to variations in terminus position (Fig. 7). In our model, the seasonal variability is likely driven by processes related to the atmospheric forcing applied (e.g. temperature and SMB variability), which in fact represents the only seasonal input used in the model. The greater-thanseasonal frequency seen in our simulations is attributed to grid resolution and missing seasonal-scale processes (e.g. ice mélange variability or seasonal ocean temperature variability) in the model. Sensitivity experiments performed on a 1 km grid did not show significant improvement with respect to ice thickness (Fig. S10) or surface speed (i.e. shape of the flow and overall magnitude; Fig. S11).
In 1990, JI had a large floating tongue (> 10 km; e.g. Motyka et al., 2011) that we are not able to simulate during the equilibrium runs. In our model (Fig. 6), the glacier starts to develop a floating tongue comparable with observations in 1999. Starting in 2000, the floating tongue is consistent in length and thickness with observations and the model is able to simulate its break-up (that occurred in late summer 2003) and the subsequent glacier acceleration. The difference between observed and modelled pre-1999 geometry results in relatively large basal melt rates for the period 1997-2003 (Fig. S9). Nevertheless, the model is able to capture the overall retreat of the terminus and the trends in the observed velocities (Figs. 2 and 3) for the period 1990-2010. Finally, the 2010-2012 observed terminus retreat (Joughin et al., 2014) is not reproduced in our simulations, likely due to inaccuracies in basal topography or misrepresentations of the atmospheric forcing and the ocean parametrization that we used. Additional sensitivity experiments showed that an increase in ocean temperature of ∼ 0.7 • C for the period 2010-2014 may trigger a retreat of the terminus that agrees better with observations (Figs. S6 and S7).
Our model reproduces two distinct flow accelerations in 1998 and 2003 that are consistent with observations. The first was generated by a retreat of the terminus and moderate thinning prior to 1998; the latter was triggered by the final breakup of the floating tongue. During this period, JI attained unprecedented velocities as high as 20 km a −1 in our simulation. Additionally, the final break-up of the floating tongue generated a reduction in buttressing that resulted in further thinning. Similar to previous studies (Nick et al., 2009;Vieli and Nick, 2011;Joughin et al. 2012), our results show that the dynamic changes observed at JI are triggered at the terminus (Figs. 7, S5, S14, and S16).
In accordance with previous studies (Thomas, 2004;Joughin et al., 2012), our findings suggest that the speeds observed today at JI are a result of thinning-induced changes due to reduction in resistive stress (buttressing) near the terminus correlated with inland steepening slopes (Figs. 6 and 7). Both model and observations suggest that JI has been losing mass at an accelerating rate and that the glacier has continued to accelerate through 2014 (Fig. 4).
The Supplement related to this article is available online at doi:10.5194/tc-10-597-2016-supplement.  -2011-IOF-301260). The development of PISM is supported by NASA grants NNX13AM16G and NNX13AK27G. We thank the editor, five anonymous reviewers for their valuable comments and suggestions to improve and clarify the manuscript, and Veit Helm for providing cryosat-2 data.