Sensitivity of the dynamics of Pine Island Glacier, West Antarctica, to climate forcing for the next 50 years

. Pine Island Glacier, a major contributor to sea level rise in West Antarctica, has been undergoing signiﬁ-cant changes over the last few decades. Here, we employ a three-dimensional, higher-order model to simulate its evolution over the next 50 yr in response to changes in its surface mass balance, the position of its calving front and ocean-induced ice shelf melting. Simulations show that the largest climatic impact on ice dynamics is the rate of ice shelf melting, which rapidly affects the glacier speed over several hundreds of kilometers upstream of the grounding line. Our simulations show that the speedup observed in the 1990s and 2000s is consistent with an increase in sub-ice-shelf melting. According to our modeling results, even if the grounding line stabilizes for a few decades, we ﬁnd that the glacier reaction can continue for several decades longer. Furthermore, Pine Island Glacier will continue to change rapidly over the coming


Introduction
Pine Island Glacier is one of the most active glaciers in Antarctica, with an ice discharge of more than 130 Gt yr −1 in 2013 Medley et al., 2014). It has experienced dramatic changes over the past decades: its velocity increased by more than 40 % between 1996 and 2007 and its grounding line retreated at a rate of about 1 km yr −1 between 1992 and 2011, which resulted in the progressive un-grounding of its ice plain (Corr et al., 2001;Mouginot et al., 2014;Rignot et al., 2014). Satellite observations reveal an average rate of mass loss multiplied by 4 between 1995 and 2006 on the main trunk (Wingham et al., 2009). The changes in ice dynamics have been attributed to the presence of warm, subsurface water in the ocean (Rignot, 1998;Payne et al., 2004), which was observed for the first time in the 1990s (Jacobs et al., 1996). The recent increase in speed was attributed to the intrusion of warm water through a widening gap in the ice shelf cavity resulting from ice shelf thinning (Jacobs et al., 2011). It was however noted that since 2009 the glacier speed at the grounding line has reached a steady value , which has been suggested to be indicative of a temporary stabilization of the glacier grounding line based on a two-dimensional model simulation (Joughin et al., 2010).
Earlier simulations with a two-dimensional model indicated a 10 % increase in velocity from a 13 % reduction in ice shelf extent, and a 70 % speed up from the collapse of the entire ice shelf (Schmeltz et al., 2002). Thomas et al. (2004b) studied the impact of grounding line migration using a flow line (1-dimensional) model. They found that the grounding line retreat of Pine Island Glacier reduced the buttressing force on the grounded part of the glacier and had a stronger effect on glacier flow than changes in ice shelf extent or thickness. They showed that grounding line perturbations were transmitted almost instantaneously over long distances inland. Their model correctly predicted that the entire ice plain would unground in the following years if ice thinning rates remain unchanged and that the ice shelf would reach a flow speed of 4 km yr −1 . Using a 2-D/3-D mixed model, Payne et al. (2004) showed that the increase in ice shelf melting would reduce basal friction at the grounding line and changes would be transmitted upstream, more than 200 km inland, on a decadal timescale, by a diffusive process. More recently, Joughin et al. (2010) used a 2-D plan-view model with constrained grounding line dynamics and ice shelf margins to conclude that the grounding line retreat would be reduced for several decades and the mass loss should remain steady. Using a volume continuity model, however, Thomas et al. (2011) found that grounding line retreat would be maintained and yield glacier speeds in excess of 10 km yr −1 within a few decades. Williams et al. (2012) concluded from a model study that high-frequency forcings (decadal to subdecadal) are transmitted by membrane stress and not by driving stress, and rapidly propagate several tens of kilometers inland. Favier et al. (2014) used a three-dimensional (3-D) full Stokes (FS) model and parameterization of the oceaninduced melting rate to study the grounding line retreat of Pine Island Glacier. They showed that the grounding line of Pine Island Glacier is likely to have started an irreversible retreat on the downward sloping bed of the main trunk and that its contribution to sea level rise could reach 100 Gt yr −1 in the next 20 yr.
Here, we use a 3-D model that includes grounding line dynamics, data assimilation for basal friction and a highresolution mesh to analyze the impact of external forcings on the ice flow dynamics of Pine Island Glacier. These external forcings are (1) surface mass balance (SMB), (2) calving front position and (3) ice-shelf melting. We discuss the impact of each external forcing on ice dynamics, i.e., on the velocity pattern over the entire basin. We compare our results with observations and conclude on the possible evolution of the glacier over the next 50 yr.

Data and methods
We initialize our numerical model to match satellite observations centered around 2008. The surface elevation is a combination of satellite radar and laser altimetry from ICESat and ERS-1 (Bamber et al., 2009). The bed topography is from ALBMAP (Le Brocq et al., 2010), derived from ground-penetrating radar measurements (Lythe and Vaughan, 2001). Ice shelf thickness is retrieved from satellite radar altimetry from ERS-1 (Griggs and Bamber, 2011) and the sea-floor bathymetry under the floating part of Pine Island Glacier is from NASA's Operation IceBridge (OIB) gravimetry . We employ surface temperature and SMB forcings (ice accumulation) from the regional atmospheric model RACMO2 (Lenaerts et al., 2012), and the geothermal heat flux is inferred from satellite magnetic data (Maule et al., 2005). The model domain corresponds to the extent of Pine Island Glacier catchment basin, which is con-strained by topography and velocity data. The grounding line position corresponds to the 2007 grounding line position derived from differential satellite synthetic-aperture radar interferometry (DInSAR) from ERS-1 and 2, RADARSAT-1 and 2 and ALOS PALSAR (Rignot et al., 2011b. We rely on melting rate reconstructions from the MITgcm using the OIB bathymetry (Schodlok et al., 2012). The melting rate is an average over a year and is kept constant throughout the simulation, so no additional melting is introduced if ice starts floating. Figure 3 shows the basal melt rate distribution under the floating ice of Pine Island Glacier as well as its distribution with depth. As bathymetric and bed data remain sparse and do not match at the grounding line, we lower the bathymetry in the first 10 km downstream of the grounding line so that hydrostatic equilibrium is consistent with the grounding line position. The correction applied to the bathymetry lowers its elevation up to 100 m over a spatially limited area. It restrains grounding line advance, which is consistent with recent observations (Rignot et al., 2011a).
We use the Ice Sheet System Model (ISSM) to perform our numerical experiments (Larour et al., 2012). A 3-D higherorder approximation (Blatter, 1995) of the full Stokes equations is applied to a 225 000 element mesh. The mesh horizontal resolution varies from 500 m near the grounding line to 10 km in the mountainous regions (see Fig. 2) and is vertically extruded in 14 non-uniform layers (thinner layers at the base). To initialize the model, the coefficient of basal friction is inferred using assimilation of InSAR-derived surface velocity data from 2008 (Rignot et al., 2011a;Mouginot et al., 2012) on grounded ice, as described in Morlighem et al. (2010) (Fig. 1a). Basal friction is assumed to follow a linear viscous law. Ice rigidity is based on the values provided in Cuffey and Paterson (2010) assuming thermal steady state on grounded ice, and is inferred using data assimilation of surface velocity on floating ice. Ice temperature and hardness are updated at each step during data assimilation of basal friction for consistency (Morlighem et al., 2010). No additional tuning, such as melting rate correction or ad hoc time-dependent friction coefficient, is applied.
The data used to initialize the model are acquired in different years, with different instruments and at resolutions that range from 300 m for observed surface velocities (Rignot et al., 2011a) to several kilometers for bedrock topography (Le Brocq et al., 2010). These data are not always consistent and lead to large ice flux divergence anomalies in ice flow simulations . We therefore relax the model for 10 yr with a fixed grounding line and use presentday forcings in order to reduce the spurious oscillations in ice thickness that exhibit large anomalies in the first years, which are caused by the uncontrolled interpolation of ice thickness data on regular grids. An alternative would be to use mass continuity to improve the bedrock topography of the grounded part of Pine Island Glacier  but it is beyond the scope of this study.  Simulations are run forward for 50 yr with time steps of three weeks to satisfy the CFL (Courant, Friedrichs et Lewy) condition (Courant et al., 1967). At each time step, the ice velocity, topography and grounding line position are updated. We use a floatation criterion for grounding line retreat: ice starts floating if it becomes thinner than the floatation thickness. Ice temperature is kept constant during the simulation and ice thickness change is computed with a mass transport equation, stabilized with the streamline upwinding finiteelement method.
We investigate the influence of external forcings through three model parameters. In the first set of experiments, we multiply the SMB (accumulation of ice) by a coefficient α varying between 0 and 3. In a second set of experiments, we simulate the retreat of the ice front position from 0 to 40 km (Fig. 1b); this is twice the distance between the 2011 ice front and the position of the rift that calved in November 2013 (Howat et al., 2012). In a final experiment, we multiply the ocean-induced melting rate pattern from the MITgcm by a coefficient β that varies between 1 and 2. A study by Jacobs et al. (2011) estimated that the ocean-induced meltwater under Pine Island ice shelf increased by 50 % between 1994 and 2009, so our multiplication factor is twice the observed increase rate over this period. In the SMB experiments, we chose a range of α such that the volume change is larger than the melting experiments. SMB, front position and melting rates are then kept constant during all the simulations. These experiments simulate changes that are twice as large as what have been recently observed.

Results
Model initialization is in good agreement with observations, with an average velocity difference of about 13 m yr −1 between modeled (Fig. 1a) and observed velocities from 2008. During the 10 yr of relaxation, the ice thickness mainly adjusts on the floating part of the glacier, with about 100 m of thickening downstream of the grounding line along the main trunk and thinning up to 150 m on the rest of the floating ice. Changes on grounded ice are more local and of smaller amplitude (Fig. 4a). Velocity is also mainly changing over the floating part of the glacier, with a speedup of 300 m yr −1 along the shear margins and slowdown of about 150 m yr −1 at the grounding line of the main trunk (Fig. 4b). A small acceleration of the main trunk and the main tributaries is also observed.  The glacier evolution during the 50 yr of the simulation under present-day conditions shows an increase in velocity over the ice shelf, from 3.7 km yr −1 in its initial state to 4.5 km yr −1 after 50 yr of evolution. This change in speed propagates several hundreds of kilometers inland: the model shows an increase in speed of 200 m yr −1 , or 20 %, 100 km upstream of the grounding line, in areas where initial speed is 1.0 km yr −1 . Changes in flow speed are detected all the way to the flanks of the glacier topographic divides; most of the glacier speeds up by 20 %. Ice thinning is equivalent to a total of 11 mm of sea level rise after 50 yr, or 78 Gt yr −1 .
The first sensitivity experiments (Fig. 5a) show that changes in SMB do not affect the ice dynamics over 50 yr simulations: ice speed changes by less that 0.1 % when SMB is multiplied by a factor of 3 compared to current values. Changes in SMB, however, lead to variations in glacier volume above floatation (Fig. 5b) equivalent to a sea level variation between −7 mm and +20 mm over 50 yr. The volume above floatation remains constant when SMB is doubled. Time series of volume change are quasi-linear because changes in volume above floatation do not involve any change in ice dynamics.
Changes in ice front position have an immediate effect on ice velocity (Fig. 5c), and make the ice front velocity vary from 3.7 to 5.4 km yr −1 . After 50 yr, the ice front velocities stabilize at the same speed, except in the case of very large retreats (>25 km, with more than half of the ice shelf being removed). Hence, changes in ice front position have only a moderate impact on long-term glacier speed. Changes in ice front position, however, do not impact the volume above floatation of the glacier (Fig. 5d). The rate of volume change is similar for all front retreat except for the most extreme retreat as all experiments lead to a similar discharge after 50 yr as noted above. The perturbations introduced here do not destabilize the glacier and do not affect its dynamic on the whole except when very large ice front retreats are introduced, as areas of the ice shelf actively buttressing the ice stream are removed in this case.
Sensitivity to basal melting under floating ice is shown in Fig. 5e and f. A doubling of the basal melting rate leads to an additional velocity increase of 800 m yr −1 on the floating part of the glacier (Fig. 5e). Acceleration is not limited to the floating part but propagates hundreds of kilometers inland in 1 to 5 yr. Increased melting also leads to a decrease in ice volume above floatation (Fig. 5f). The time series of volume above floatation show an increased contribution with time as the glacier is accelerating in response to changes in basal melting rate. Multiplying the basal melting by a factor of 2 leads to an additional ice volume loss of 4 mm of sea level equivalent.  In all the above scenarios, variations in the grounding line position are small and grounding line positions remain in areas with a fine mesh resolution. Figure 6a shows the grounding line position at the beginning of the simulations and at the end of the melting scenarios for β = 1 and β = 2. The grounding line retreats by no more than 10 km during the 50 yr simulations and is very limited over the main trunk. Using a finer mesh resolution did not change the results, proving that grounding line retreat is not limited by a too-coarse mesh resolution (see appendix and between the basal melting rates that we are applying (based on high-resolution MITgcm simulations), which are significantly lower that the ones used in Favier et al. (2014).

Discussion
Changes in both SMB and basal melting affect Pine Island Glacier's volume, but basal melting under floating ice is the only modeled forcing that affects the glacier's dynamics on the timescales under study here. Increased basal melting causes thinning of floating ice that leads both to a reduction in buttressing from the ice shelf and a grounding line retreat. Experiments focusing on ice front retreat also confirm that limited ice front retreat over an unconfined part of the ice shelf, due to calving events for example, have no long-term effect on the glacier's dynamics.
In their study of Pine Island Glacier, Favier et al. (2014) show that the grounding line of Pine Island Glacier is engaged on an unstable 40 km retreat and that the glacier is controlled by marine ice sheet instability. Their results also show that limited ice front retreats do not affect grounding line dynamics, while changes in basal melting rates under floating ice strongly impact grounding line motion. In their control experiment, the basal melting is parameterized to match certain recent observations (Dutrieux et al., 2013). A small reduction in grounded ice area is observed in this case, which is similar to the results reported here. Their simulations show that if basal melting increases and extends to a larger portion of the ice shelf, the grounding line starts an unstable retreat along the 40 km retrograde slope. In our simulations, even in the case of doubled melting rate, the grounding line position does not retreat more than 10 km. This is probably caused by the different patterns of melting rates: basal melting rates in  (Fig. 10). This confirms that the applied basal melting rate is critical for simulations of ice stream dynamics and grounding line retreat.
Our simulations reveal that even if increased basal melting causes limited grounding line retreat, it reduces the buttressing from the ice shelf as the ice is thinning, which leads to a speedup of Pine Island Glacier. A change in basal melting not only affects ice velocity in the floating part of the glacier: acceleration propagates inland, and reaches the flanks of the ice divide, as predicted by Williams et al. (2012) for decadal forcings. Our simulated accelerations are difficult to compare with previous results, as acceleration is rarely provided. However, they propagate further inland than in most prior studies; we obtain a velocity increase of about 100 m yr −1 up to 200 km upstream of the grounding line, through a combination of transmission of membrane stress and driving stress or diffusive processes. In previous studies, a similar speedup is confined to the first 70 km upstream of the grounding line in Payne et al. (2004) and to the first 120 km in Thomas et al. (2004a). In Joughin et al. (2010), the thinning rate propagates inland and reaches 1 m yr −1 100 km upstream of the grounding line after a decade.
Comparison of the first 15 yr of simulation with satellite observations of previous acceleration of Pine Island Glacier in the 1990s and 2000s   (Fig. 7d-f) provides qualitative estimates, as model and observations are from different years. It shows that the patterns of modeled acceleration due to enhanced sub-ice-shelf melting (α = 1.5) are in agreement with the observed glacier acceleration, with similar patterns after 10 and 15 yr. Modeled velocities differ from observations along the side margins of the ice shelf: in this region, ice accelerates more in the model than in the observations (1000 m yr −1 in the first 10 yr of simulation vs. 800 m yr −1 between the 1996 and 2006 observations). Our simulation shows an acceleration of the main trunk and most of its tributaries (Fig. 7b, e)  plain after 15 yr of simulations vs. 900 m yr −1 in the observations, Fig. 7b, e), suggesting that our results underestimate the actual speedup of Pine Island Glacier. This acceleration is difficult to compare to that inferred in prior studies, which mainly focused on the glacier centerline.
In an additional experiment (Fig. 8a) we increase sub-iceshelf melting for a limited time. In this simulation, basal melting from the MITgcm is multiplied by 1.5 for the first 5 or 15 yr and then switched back to its initial value. The ice shelf velocity increases and reaches ∼ 4.9 km yr −1 in both cases after 50 yr compared with 5.0 km yr −1 when increased basal melting (also multiplied by 1.5) is kept constant over the entire 50 yr simulation. In the control run where basal melting from MITgcm is directly used, the ice shelf velocity after 50 yr is 4.5 km yr −1 . This indicates that a temporary increase in basal melting rates has a long-term impact on ice dynamics and that Pine Island Glacier will not slow down and stabilize, even if ocean conditions were to return to what they were a few decades ago. This conclusion is consistent with the marine instability hypothesis that states that glaciers on downward-sloping bed inland are intrinsically unstable (Weertman, 1974;Schoof, 2007) and with recent studies of Pine Island Glacier (Favier et al., 2014). Change in glacier volume above floatation after 50 yr is almost identical if the basal melting rate is increased for 5 or 15 yr (Fig. 8b), and is about 0.4 mm of sea level equivalent lower than if increased melting is kept constant for 50 yr.
Our study provides estimates of climate sensitivity of Pine Island Glacier based on a 3-D higher-order formulation, with a high resolution in the grounding line area. No melting rate or surface accumulation correction is introduced to start with a model in steady-state condition (Joughin et al., 2010;Cornford et al., 2013), and no additional parameterization is needed to include buttressing (reduction of basal friction), contrary to most of the studies performed with flow band models. We have shown in another study that errors in  ice rigidity and basal friction do not affect the results significantly for these short-term simulations . However, our model has some limitations, such as a fixed ice front that can only be changed manually, no rheological weakening of the floating ice and a grounding line that is not based on contact mechanics, which would be too computationally intensive for this kind of sensitivity experiments.
In all our simulations, the pattern of basal melting is kept constant with time. Additional experiments have been performed to show that introducing moderate melting rates under ungrounding ice does not significantly affect our results. Figure 9 shows the ice velocity and ice volume above floatation if we apply either the average basal melting rate under floating ice or the highest value of melting rate in the areas that unground. Results show that applying the average basal melting rate value has a small impact on the simulations, while extreme melting leads to variations comparable to a doubling of the basal melting rate over the floating area. This result is not surprising, as the grounding line retreat remains limited in all our experiments; however, the results could be different in the case of a large grounding line retreat. This is comparable to results of Favier et al. (2014) based on Elmer/Ice; they run two experiments with similar basal melting parameterization. In the first one, basal melting is turned on as grounded ice starts floating, while in the second one, basal melting is limited to the initial floating part of the glacier. Results of these two experiments are very similar in terms of both grounded area and sea level rise.
Melting rates are kept constant throughout the simulations, while we know that changes in ice shelf cavity will affect their amplitude and spatial distribution (Schodlok et al., 2012). We choose not to change the pattern of basal melting, as we do not know how changes in ice shelf cavity will impact oceanic circulation and basal melting rates, and our results are therefore likely to be conservative estimates of changes. Melting rate parameterizations (Pollard and De-Conto, 2009;Little et al., 2009) provide a first estimate but do not include specifics for each glacier. Results from the MITgcm highlight that no simple parameterization of basal melting rate based on ice shelf depth, for example, can be derived (Fig. 3). Recent observations in the bay adjacent to Pine Island Glacier also report the large temporal variability of ocean heat and ocean-induced melting in this area (Dutrieux et al., 2014). Our results show that precise estimates of basal melting under floating ice are required and essential for constraining the evolution of glacier dynamics, as a modest increase of 10 % in basal melting rates impacts ice sheet dynamics. To achieve this goal, however, progress is necessary in the modeling of ice-ocean interactions beneath the ice shelves with coupled ice-sheet-ocean-atmosphere models (Schodlok et al., 2012). Finally, our simulations suggest that the acceleration of Pine Island Glacier will continue to propagate inland and its mass loss will continue for decades to come, even if the oceanic conditions return to their state prior to the 1990s and the grounding line position remains stable for a few decades. Similarly, if more ocean heat reaches the grounding line area, the mass loss will continue to increase for decades to come, with no sign of stabilization.

Conclusions
Our study shows that Pine Island Glacier is highly sensitive to basal melting under its floating extension; this parameter is the most important climate factor tested in our simulations and controls most of the dynamics of this glacier, even if grounding line retreat remains limited. The effect of changes in sub-ice-shelf melting are not limited to the floating tongue and grounding line area but are rapidly transmitted hundreds of kilometers inland. Increase in sub-ice-shelf melting for only 5 yr destabilizes the glacier for several decades and has a long-term impact on its dynamics. A comparison of our model results with satellite observations of the 1990s and 2000s shows the good qualitative agreement between modeled and observed accelerations and suggests that the glacier speedup is consistent with increased basal melting under the ice shelf, although simulations and observations are separated by a decade. Coupling of the ice sheet with ocean circulation models is therefore desired for future studies to conduct more accurate simulations, as the glacier is controlled by ocean-induced melting rates. Overall, Pine Island Glacier is likely to keep accelerating over the coming decades, even if ocean circulation and melting remain unchanged and the grounding line temporarily stabilizes over a local feature in the bedrock topography.