Glacial isostatic adjustment in response to changing Late Holocene behaviour of ice streams on the Siple Coast, West Antarctica

The Siple Coast region of Antarctica contains a number of fast-ﬂowing ice streams, which control the dynamics and mass balance of the region. These ice streams are known to undergo stagnation and reactivation cycles, which lead to ice thickness changes that may be sufﬁcient to excite a viscous solid Earth response (glacial isostatic adjustment; GIA). This study aims to quantify Siple Coast ice thickness changes during the last 2000 yr in order to determine the degree to which they might contribute to GIA and associated present-day bedrock uplift rates. This is important because accurate modelling of GIA is necessary to determine the rate of present-day ice-mass change from satellite gravimetry. Recently-published reconstructions of ice-stream variability were used to create a suite of kinematic models for the stagnation-related thickening of Kamb Ice Stream since ∼ 1850 AD, and a GIA model was used to predict present-day deformation rates in response to this thickening. A number of longer-term loading scenarios, which include the stagnation and reactivation of ice streams across the Siple Coast over the past 2000 yr, were also constructed, and used to investigate the longer term GIA signal in the region. Uplift rates for each of the ice loading histories, based on a range of earth models, were compared with regional GPS-observed uplift rates and an empirical GIA estimate. We estimate Kamb Ice Stream to have thickened by 70–130 m since stagnation ∼ 165 years ago. Modelled present-day vertical motion in response to this load increase peaks at − 17 mm yr –1 (i.e. 17 mm yr –1 subsidence) for the weakest earth models tested here. Comparison of the solid Earth response to ice load changes throughout the last glacial cycle, including ice stream stagnation and reactivation across the Siple Coast during the last 2000 yr, with an empirical GIA estimate suggests that the upper mantle viscosity of the region is greater than 1 × 10 20 Pa s. When upper mantle viscosity values of 1 × 10 20 Pa s or smaller are combined with our suite of ice-load scenarios we predict uplift rates across Siple Coast that are at least 4 mm yr –1 smaller than those predicted by the empirical GIA estimate. GPS data are unable to further constrain model parameters due to the distance of the GPS sites from the study area. Our results demonstrate that Late Holocene ice load changes related to the stagnation and reactivation of ice streams on the Siple Coast may play a dominant role in deﬁning the present-day uplift signal. However, both the detailed Earth structure and deglacial history of the region need to be better constrained in order to reduce uncertainties associated with the GIA signal of this region.

. Study region with rock outcrops in brown. Profile locations 1-5 from Retzlaff & Bentley (1993) shown as red dots. Masks indicate the previously fast-flowing Kamb Ice Stream (KIS, blue outline), the extent of the KIS ice loading calculation (orange outline) and Whillans (WIS) and MacAyeal (MacIS) ice streams (green shaded areas). GPS locations are shown as purple circles. The Duckfoot region is shaded in grey.
in this region over the past ∼1000 yr (e.g. Hulbe & Fahnestock 2007;Catania et al. 2012), with ice streams displaying asynchronous century-scale stagnation and reactivation. The focus of this study is to explore how such century-scale changes might have affected the glacial isostatic adjustment (GIA) signal across the Siple Coast, and hence regional estimates of present-day ice mass change derived from Gravity Recovery and Climate Experiment (GRACE) data. Isolation of that component of the GRACE signal which is due to present-day ice mass change relies on an accurate estimate of the signal due to past ice mass change, that is GIA.
Ice streams are fast flowing due to the presence of water at their bed and deformable sediments below. It is thought that stagnation occurs due to an increase in basal resistance related to strengthening of the underlying sediments, caused by basal freezing or a decrease in subglacial water pressure (Beem et al. 2014). Kamb Ice Stream (KIS, formerly 'Ice Stream C', see Fig. 1) stagnated around 1850 AD (Retzlaff & Bentley 1993) or perhaps even earlier (Shabtaie & Bentley 1987). Altimeter measurements of surface elevation change for the periods 2003-2007 and 2011-2014 show that the trunk of KIS is currently thickening by up to 0.6 m yr -1 (Pritchard et al. 2009;Helm et al. 2014) as ice continues to flow into the drainage basin from upstream but no longer flows out (Rignot et al. 2011). GPS observations of neighbouring Whillans Ice Stream (WIS) show it is decelerating, and given the long-term (decadal) average rate of ice flow it has been suggested it may stagnate in the next ∼50-150 yr (Joughin et al. 2005), or possibly sooner if deceleration continues to increase (Beem et al. 2014).
It is not known how changes in ice thickness due to the stagnation and reactivation of ice streams on the Siple Coast could affect GIA in this region. Existing forward models of GIA, which generally do not consider ice thickness changes during the last few millennia, tend to predict high uplift rates along the Siple Coast [e.g. 10 mm yr -1 , Whitehouse et al. (2012b); 11 mm yr -1 , Argus et al. (2014)], but there are few independent geological data to constrain GIA models in this region (Whitehouse et al. 2012a) and bedrock GPS sites are remote due to a lack of rock outcrops (Fig. 1).
A recent study by King et al. (2012) examined GRACE time series over the Ross Ice Shelf and suggested that the GIA signal in this region should be close to zero to avoid implausible estimates of substantial mass loss over the ice shelf. This argument is based on an understanding that GRACE-derived mass changes over the floating Ross Ice Shelf should be free of ice-mass change signal (once signal leakage is considered) due to hydrostatic equilibrium; as such the GRACE signal should be dominated by GIA and small ocean mass changes. However, King et al. (2012) found that leakagecorrected GRACE mass change over the area-averaged Ross Ice Shelf already had near-zero trend, implying that the GIA signal in this region was negligible. Because the gravitational effect of GIA is long-wavelength it was inferred that present-day Siple Coast GIA was likely also small. King et al. (2012) consequently set the GIA correction to zero for two Siple Coast ice drainage basins, resulting in changes of ∼±6-8 Gt yr -1 in estimated ice mass trends. This result is in stark disagreement with existing GIA model predictions, and highlights the need for improved understanding of the retreat history of this region.
Ice stream stagnation-reactivation will cause an increase/decrease in ice thickness on century timescales, and, given a sufficient amount of thickening/thinning, there may be a substantive GIA-related response depending on the local properties of the Earth. Most GIA models neglect ice load changes across Antarctica during the last 1000 or 2000 yr (e.g. Whitehouse et al. 2012b;Ivins et al. 2013;Argus et al. 2014) and hence will not account for ice thickness changes due to recent ice stream re-organisation along the Siple Coast. The exceptions to this are GIA models that make use of a transient ice-sheet reconstruction, derived by coupling a numerical ice-sheet model to a sea level solver (e.g. Gomez et al. 2013;de Boer et al. 2014). Given the difficulty of reproducing ice stream stagnation and reactivation within a dynamic model (e.g. Kirchner et al. 2011), and the sensitivity of the modelled Siple Coast grounding line position to ice-ocean feedbacks (Gomez et al. 2013), it is unlikely that these coupled models accurately capture the details of Siple Coast ice stream dynamics during the last thousand years. Moreover, these models used an upper mantle viscosity that may be too high for this region (van der Wal et al. 2015), which would lead to an underestimate of the decadal response to ice thickness changes. Despite these limitations, such models have the potential to provide an important insight into some of the processes that control ice stream dynamics along the Siple Coast, where the low-gradient bathymetry means that long-wavelength changes in the height of the solid Earth or the shape of the geoid may have a profound effect on grounding line and ice shelf dynamics (Greischar & Bentley 1980;Gomez et al. 2010).
In this study, the potential perturbation to the present-day GIA signal due to the stagnation and reactivation of ice streams, and the consequent thickening and thinning of the ice sheet over the last 2000 yr, is investigated. A specific goal of the study is to test the hypothesis that GIA-related subsidence in response to recent ice stream re-organization may be sufficient to reduce GIA modelpredicted uplift rates along the Siple Coast to negligible levels, as suggested by King et al. (2012). Given the absence of rock outcrops within this region we compare our revised GIA modelpredicted uplift rates with independent empirical estimates based on a combination of GRACE, satellite altimetry and a firn densification model (Gunter et al. 2014).
The following section describes the methods used to (i) reconstruct the ice loading history of the Siple Coast (Section 2.1) and (ii) determine the present-day solid Earth response to these load changes (Section 2.2). The data sets that were used to constrain the modelling are described in Section 2.3. The results of the ice sheet reconstructions and the GIA modelling are presented in Section 3, and these are discussed further in Section 4.

Ice loading history reconstructions
Due to the limitations of ice-sheet models in representing ice streams in this region (e.g. Pollard & DeConto 2012;Whitehouse et al. 2012a), an ice-sheet model was not used to generate the recent Siple Coast loading history. Instead, a simpler semi-empirical approach was used, based only on available observations relating to ice input and output across the Siple Coast. Sufficient observations exist to permit the detailed reconstruction of surface load changes associated with the most recent stagnation of KIS ∼165 years ago (see Section 2.1.1). Prior to this, the existence of stagnation-reactivation events for several of the Siple Coast ice streams has been inferred from flow features on the Ross Ice Shelf (Hulbe & Fahnestock 2007;Catania et al. 2012), but constraints on the timing and magnitude of these events are limited. Hypothesized surface load changes associated with these earlier stagnation-reactivation events are described in Section 2.1.2., and a summary of the ice load histories used to drive the GIA model is listed in Section 2.1.3.
It is important to note that our intention is not to reconstruct an accurate ice sheet history consistent with ice dynamics, but to consider a range of plausible ice thickness change scenarios, based on observations and models available in the literature, for the purpose of investigating the solid Earth response to surface loading. In order to set up our experiments we reconstruct temporal variations in ice thickness across the region of stagnation by considering ice mass inputs and outputs to the system. We have attempted to account for uncertainties associated with the timing and magnitude of past ice flow changes when creating our ice thickness reconstructions.
The rate of ice input into the stagnation region will depend on the evolution of the regional ice velocity field as well as local surface mass balance rates. Past ice velocities are poorly known (and modelling the evolution of the ice sheet is outside the scope of this study), so we address this uncertainty by using a range of estimates for the pre-stagnation velocity of KIS. We assume that ice flow occurred via basal sliding, that is that velocities were constant throughout the depth of the ice stream, and that upstream velocities were constant throughout the period of stagnation.
The reality will have been far more complex (e.g. Price et al. 2001;Conway et al. 2002). Feedbacks between ice dynamics, subglacial hydrology, basal melt rates and till strength are hypothesized to have influenced ice stream stagnation and flow-switching via a variety of processes (e.g. Tulaczyk et al. 2000;Hulbe & Fahnestock 2004;Joughin et al. 2004;van der Wel et al. 2013), including water piracy (Anandakrishnan & Alley 1997) and the lateral migration of shear zones (Jacobson & Raymond 1998). Whatever processes were responsible for centennial-scale changes across the Siple Coast, they will have resulted in a complex pattern of spatial and temporal variations in thickening and thinning of the ice sheet. Our experiments are a first-order attempt at reconstructing representative ice thickness changes across the region. Alternative reconstructions should be explored as our understanding of the dynamics of this system improves. In particular, there are many uncertainties associated with the timing and magnitude of the oldest stagnation and reactivation events that we consider (see Section 2.1.2). Any assumptions we make when constructing our ice histories, or when using data to constrain our calculations, are highlighted in the appropriate section and discussed further in Section 4.

KIS stagnation
The build-up of ice associated with the stagnation of KIS is due to the imbalance of ice input and ice output. It is assumed that before the commencement of stagnation the ice stream was in balance that is any input was balanced by flow. As the ice stream stagnated, it thickened due to a combination of continued ice flux from upstream, local accumulation, and a lack of flow downstream. Therefore, in order to calculate the build-up of ice, the velocity of the ice stream before stagnation is required, as well as accumulation rates since the initiation of stagnation (Section 2.3.1). Present-day rates of ice elevation change also provide an estimate of post-stagnation thickening rates (Section 2.3.2).
The time since stagnation also affects the amount of ice build-up, with earlier estimates for the initiation of stagnation resulting in greater estimates for the total amount of ice thickening. There are uncertainties associated with the timing of stagnation and the speed at which KIS was flowing prior to stagnation (see Sections 2.1.1.1 and 2.1.1.2) so, to take account of these, several stagnation histories were constructed to explore the parameter space of the variables.

Pre-stagnation velocity on KIS.
The stagnation of KIS is reported to have occurred approximately 165 years ago (Retzlaff & Bentley 1993) and prior to this KIS is thought to have been fast flowing, although the pre-stagnation velocity is uncertain. Ng & Conway (2004) estimate that KIS flowed at an average of 340 m yr -1 over the 740 yr before stagnation, while a modelling-based study by Hulbe & Fahnestock (2007) used values of 300 and 500 m yr -1 for lower and upper bound velocities, respectively, for flow ∼200 yr before stagnation. It is thought that KIS stagnation was initiated by the rapid stagnation of an area known as 'Duckfoot' (see Fig. 1) around 360 years ago (Catania et al. 2012), which eventually led to the stagnation of the main trunk of the ice stream around 200 yr later. Following the shutdown of the Duckfoot region, ice velocity at the margins of KIS is thought to have slowed to ∼200 m yr -1 (Catania et al. 2006), although the authors note that this is likely to be an underestimate of the centreline velocity.
In order to calculate the change in ice thickness during the stagnation of KIS it is necessary to know the ice velocity immediately prior to stagnation of the main trunk, that is after stagnation of the Duckfoot region and the consequent narrowing and slowing of the main trunk. Uncertainty in the pre-stagnation velocity is accounted for by using three different values that capture the full range of plausible velocities: 100, 200 and 300 m yr -1 . 300 m yr -1 was used as a lower bound by Hulbe & Fahnestock (2007), but that modelling study referred to the period prior to stagnation of the Duckfoot region. Here, this value is used as an upper bound on the pre-stagnation velocity for KIS immediately prior to shutdown of the main trunk.
2.1.1.2 Timing and rate of stagnation on KIS. The timing of KIS stagnation is uncertain with studies suggesting stagnation times 100-150 yr apart (e.g. Shabtaie & Bentley 1987;Retzlaff & Bentley 1993). We base our reconstructions on the estimated age of buried crevasses reported by Retzlaff & Bentley (1993), who were able to determine the timing of stagnation at five profiles across KIS using ice-penetrating radar to identify buried crevasses (see Fig. 1). These authors also provide uncertainty bounds for the timing of stagnation at each profile, with all ages being given in terms of 'years ago.' Here, timings are reported in calendar years, and 'years ago' is taken to refer to the number of years prior to 1988 [the year in which the radar data discussed in Retzlaff & Bentley (1993) were collected]. The timing of stagnation at each profile is given in Table 1. Shabtaie & Bentley (1987) suggest an earlier stagnation time than Retzlaff & Bentley (1993), by around 100 yr (equivalent to a date of 1735 AD at the location of profile 3). However, since they estimate stagnation time at one location only a stagnation history cannot be built up from their data alone. Instead we choose to use the later timings given by Retzlaff & Bentley (1993), which provide a lower bound on net ice build-up and hence the magnitude of the solid Earth response. We discuss the implications of this uncertainty in stagnation timing in Section 4.
It is thought that the stagnation of KIS occurred in a wave moving upstream from the grounding line (Retzlaff & Bentley 1993). Profiles 1-3, lying closest to the grounding line, have the same Table 1. Timing of stagnation of KIS, based on information in Retzlaff & Bentley (1993).

Profile
Timing of stagnation  Retzlaff and Bentley (1993); we assume KIS took 5 yr to stagnate from the grounding line to profile 1, see Table 2. b LB, BE and UB timings explore the range of times at which profiles 1-3 are thought to have stagnated (1828-1888). c Uncertainty not given. d The rate at which stagnation proceeds is defined to be the same in all experiments (final column), and hence the total time taken for stagnation between the grounding line and profile 5 is the same in each case. These constraints preclude us from exploring the full suite of hypothesized stagnation times at profile 5 (1888-1958). estimated age of stagnation (Table 1), so we used uncertainty bounds of ±30 yr (Retzlaff & Bentley 1993) to create three ice loading histories with different timings. The total time taken for stagnation to occur between the grounding line and profile 5 is the same in all loading histories, but in our end member cases stagnation commences earlier (upper bound) and later (lower bound) than the best estimate. Since ice continues to accumulate following stagnation (see Section 2.1.1.4), these end member scenarios allow more or less ice to build up between the end of stagnation and the present day, respectively. In the upstream reaches of the ice stream, observations of presentday ice velocity show that tributaries are still flowing into the trunk of KIS at approximately 50 m yr -1 (Rignot et al. 2011). In fact, ice is still flowing up to the location of profile 4, so we treated the ice thickness change calculation between profiles 4 and 5 differently to the downstream sections (see Section 2.1.1.4).
The speed at which stagnation propagates upstream is determined by the distance between the profiles and the timings in Table 1. Since Retzlaff & Bentley (1993) report that stagnation occurred rapidly, and the first three profiles have the same stagnation age, we kept the rate of propagation the same for all variations of the ice loading history. This rate decreases upstream, as shown in Table 1.
It is likely that at a given location the ice stream velocity decreased over a number of years prior to complete stagnation. However, for ease of computation we assumed that the stagnation of each 5 km grid cell (see Section 2.1.1.4) occurred in 1 yr. After a section of the ice stream has stagnated, we assumed that its velocity is zero, so that no ice flows into or out of this location. The set of reconstructed loads should capture the full range of loading scenarios, as we have adopted conservative bounds in the calculations.

Spatial extent of KIS stagnation.
In order to isolate the response to KIS stagnation, our GIA model limits ice thickness changes during the last ∼165 yr to the KIS drainage basin, with thickness change values set to zero everywhere else. By limiting the study region to the KIS basin for this period, any GIA-related deformation due to thickening or thinning of neighbouring ice streams will be neglected, but the impact on our results is likely to be minimal as KIS is the only ice stream observed to experience significant thickening over the past ∼200 yr (Catania et al. 2012). Ice streams to the north of KIS appear to be in near-equilibrium, with presentday elevation change rates less than ±0.1 m yr -1 , which can be attributed to variations in surface mass balance (McMillan et al. 2014). Neighbouring WIS is currently thinning in its upstream areas and thickening near the grounding line (Pritchard et al. 2009), but these present-day changes have only begun recently (Joughin et al. 2005) and hence were not included in this component of the study. Changes to WIS are included our reconstruction of the longer term (∼2000 yr) ice loading history of the region in Section 2.1.2.
The area over which the KIS ice loading history was calculated is shown by the orange outline in Fig. 1. This outline was created by expanding the drainage basin from Rignot et al. (2008) to include the tributaries still flowing into KIS and the region of present-day thickening shown by Pritchard et al. (2009). The drainage basins derived by Zwally et al. (2012) were also considered, but not used as the drainage basin for KIS also includes a part of WIS in that data set. The area of fast ice stream flow prior to stagnation was defined based on fig. 3 of Catania et al. (2012), and is shown in blue on Fig. 1. In defining this region it is assumed that fast-flowing ice extended upstream to the location of profile 5, which coincides with the extent of fast flow in neighbouring ice streams (Rignot et al. 2011).

Net ice thickness change.
We reconstructed KIS ice thickness change on a 5 × 5 km grid for two periods: during stagnation, and from the end of stagnation to present. During the first period, ice thickness changes were calculated on a yearly basis. It is assumed that upstream of stagnation the ice sheet is in balance, that is, accumulation is balanced by flow so that there is no change in ice thickness. At the stagnating grid cells, net ice-thickness change is given by the sum of the influx of ice plus accumulation. Grid cells downstream of stagnation, where the ice flux is zero, are assumed to be out of balance, with accumulation causing ongoing ice-sheet thickening.
When calculating the influx of ice into a cell it is assumed that flow occurs by basal sliding and hence ice velocity is constant throughout the full depth of the ice stream. Influx for a given cell within the ice stream area is therefore: Ice thickness is taken from Bedmap2 (Fretwell et al. 2013). It is assumed that any differences between the recent (Bedmap2) and pre-stagnation ice thicknesses will have a negligible effect on our results, and this is discussed further in Section 4.2. Furthermore, we assumed that ice-thickening would not have been confined to the area of previously fast-flowing ice, and in reality adjacent regions would have accommodated some of the influx. This was accounted for by assuming that in any given year 80 per cent of the total influx contributed to thickening grid cells in the previously fastflowing region (blue outline, Fig. 1), while 20 per cent contributed to thickening neighbouring grid cells of the drainage basin (orange outline, Fig. 1). Influx was divided by the area of the grid cells to obtain the thickening in each region and this field was smoothed using a Gaussian filter (with 17.5 km half width). In addition to influx, accumulation also contributes to net ice-thickness change for stagnating or stagnated (downstream) cells as there is little or no outflow from these regions. As mentioned previously, ice velocity is assumed to be zero downstream of the propagating wave of stagnation. However, between profiles 4 and 5, ice currently still flows at ∼50 m yr -1 (Rignot et al. 2011). Ice velocity in this region was therefore assumed to linearly decrease from its pre-stagnation value to the presentday value (50 m yr -1 ) over the period indicated in Tables 1 and 2. Upstream of profile 5 it is assumed that the ice sheet is in balance during stagnation (no net thickness change). Full details of the calculations used to reconstruct ice thickness changes during stagnation are given in Table 2.
In order to determine the full present-day GIA response to KIS stagnation, ice thickness changes between the completion of the stagnation (see Table 1) and the present day must also be accounted for. In the upper reaches of KIS, post-stagnation ice thickness change due to the ongoing influx of ice from upstream is deduced from ICESat elevation change rates where these are greater than 0.3 m yr -1 (as shown by the dashed contour in Fig. 2). Following McMillan et al. (2014), we assume any elevation changes <0.3 m yr -1 in the ICESat data are purely due to short-term fluctuations in surface mass balance (SMB), a signal that should not be extrapolated over the full period since the completion of stagnation. Therefore, outside the 0.3 m yr -1 contour we assume thickening is due to accumulation only and use the long term average SMB rates from RACMO2.1/ANT .  Table 2. Details of the calculations used to reconstruct ice thickness changes during stagnation. 'Profiles' refer to locations described by Retzlaff & Bentley (1993).

Location Description
Grounding line to Profile 1 Profile 1 is the furthest downstream and is around 25 km from the current grounding line. It is assumed that stagnation between the grounding line and profile 1 occurred over 5 yr at a rate of 5 km yr -1 . The initial time of stagnation for a range of scenarios is given in Table 1.
Profiles 1-4 Between profiles 1 and 4 ice-sheet thickness change is calculated on 1 yr time steps. As stagnation propagates upstream, the number of grid cells 'stagnating' at a given time depends on how far upstream stagnation has reached and the speed of stagnation as detailed in Table 1. Further details are given below.
Profiles 1-2 The ice stream stagnated at a rate of 15 km yr -1 , so 3 grid cells in the along-stream direction stagnate in each 1 yr time step, based on 5 km grid cells.
Profiles 2-3 The ice stream stagnated at a rate of 10 km yr -1 -2 grid cells per time step.
Profiles 3-4 The ice stream stagnated at a rate of 5 km yr -1 -1 grid cell per time step.
Profiles 4-5 The section between profiles 4 and 5 is treated differently because ice is still flowing here. In this section, flow of the ice stream is slowed down rather than stagnated completely. The velocity is reduced from the pre-stagnation values (i.e. 100, 200 and 300 m yr -1 ) to the present-day velocity (50 m yr -1 ) linearly over the length of time this section takes to stagnate (28 time steps, see Table 1). A major assumption with this approach is that the pattern of thickening observed by ICESat between 2003 and 2007 has remained the same since the ice stream stopped stagnating. It should also be noted that ICESat observations relate to surface elevation change rather than ice thickness change, that is they may contain a component due to bedrock motion. This discrepancy is not accounted for here because bedrock rates are likely to be at least an order of magnitude smaller than the ICESat cut-off value of 0.3 m yr -1 Whitehouse et al. 2012b), and therefore the error incurred by using uncorrected ICESat rates to calculate   Table S1) ice thickness change will have a negligible effect on GIA model output.
In total, nine ice loading histories for KIS stagnation have been constructed by combining the three timing scenarios detailed in Table 1-lower bound (LB), best estimate (BE), and upper bound (UB)-with each of the three pre-stagnation velocities discussed in Section 2.1.1.1 (100, 200 and 300 m yr -1 ). The choice of combination is used to define the name of each model, for example LB_100 refers to the lower bound timing scenario combined with a pre-stagnation velocity of 100 m yr -1 .

Regional loading history
In addition to the standalone KIS stagnation models described in Section 2.1.1, a loading scenario has been modelled which includes a 2000-yr history of stagnation and reactivation of KIS, WIS and MacAyeal Ice Stream (MacIS). Recent load changes to these ice streams are combined with the W12 ice loading history (Whitehouse et al. 2012a) in order to create a model of ice load changes throughout the last glacial cycle, noting that the W12 model adopts the ICE-5 G deglacial model outside Antarctica (Peltier 2004), and does not include any ice load changes during the last 2000 yr.
The Siple Coast ice stream history summarised by Catania et al. (2012) was used to construct a simple loading history which accounts for stagnation of the WIS and MacAyeal Ice Stream (MacIS) 850 yr BP (before present), and their reactivation 400 yr later. As detailed data are not available regarding the timing of the stagnation of these ice streams, a simple model was constructed using a uniform amount of ice thickening (200 m) over a small portion of the ice stream, as shown by the green shaded regions in Fig. 1. These areas represent the currently fast-flowing sections of WIS and MacIS where maximum ice build-up might be expected following stagnation. Ice loading and unloading was applied instantaneously at the start and end of the 400-yr period, respectively (see Table 3).
The Siple Coast ice streams have been predicted to display cyclic behaviour (Robel et al. 2013). Therefore, in addition to the load changes on WIS and MacIS described above, an earlier stagnationreactivation cycle of KIS is also included. Ice loading and unloading associated with this event are applied instantaneously at 1800 and 1000 yr BP, respectively, and the magnitude of ice thickening and thinning is taken to be equivalent to the net ice thickness change in the BE_200 (KIS-only) ice model (see Section 2.1.1.4).
The final loading event in our regional reconstruction is the progressive build-up of ice associated with the most recent stagnation of KIS. Ice load changes are applied at a series of discrete times, as per model BE_200. This model is chosen from the suite of nine stagnation models because it adopts the most likely pre-stagnation velocity and stagnation duration.
All of these post-2 ka BP ice load changes are added onto the end of the W12 deglaciation model (see Table 3 for detailed timings). The W12 deglacial history is used rather than the W12a deglacial history (Whitehouse et al. 2012b) because the latter model includes the addition of significant amounts of ice to the Antarctic Peninsula during the last 1000 yr, which may affect the GIA signal in the Siple Coast when combined with a weak earth model. It is likely that ice thickness changes along the Antarctic Peninsula during the last 1000 yr were in fact an order of magnitude smaller than those included in the W12a model, and hence are unlikely to have influenced Siple Coast deformation (Nield et al. 2012;Wolstencroft et al. 2015).

Ice load scenarios for input to the GIA model
In total, three ice load scenarios have been used to drive the GIA model: (1) The KIS stagnation scenarios described in Section 2.1.1 are used in isolation to investigate the solid Earth response to this most recent loading event. These models are referred to as KIS-only models.
(2) The regional scenario described in Section 2.1.2, which combines the Antarctic-wide W12 ice load model with a detailed history of Siple Coast ice loading during the last 2000 yr, is used to investigate the solid Earth response to recent ice load changes in the context of longer-term deglaciation. This is referred to as the W12+Siple model, and will be compared with previously-published results associated with the W12 ice load model (Whitehouse et al. 2012b).
(3) An additional ice loading history is constructed in which grounding line retreat in the Ross Sea component of the W12 model is accelerated such that the present-day ice extent and thickness is reached by 5 ka BP rather than 2 ka BP. This present-day ice configuration is then held constant between 5 and 2 ka BP before the regional loading history for the Siple Coast is invoked during the last 2000 yr.
Scenario 3 explores the hypothesis that the bulk of Ross Sea grounding line retreat occurred during the early-to-mid Holocene. The majority of evidence relating to grounding line retreat in this region comes from the western coast of the Ross Sea, adjacent to the Transantarctic Mountains, and it suggests that retreat persisted into the Late Holocene (e.g. Hall et al. 2013;Anderson et al. 2014). However, the complex topography along this coast provides pinning points that will have delayed grounding line retreat in this area relative to the rest of the Ross Sea. We therefore test a scenario in which grounding line retreat across the whole of the Ross Sea is at University of Durham on March 31, 2016 http://gji.oxfordjournals.org/ Downloaded from complete by 5 ka BP. This early retreat model is referred to as the W12_5k+Siple model. Details of our reconstructed ice thickness change across KIS, MacIS and WIS are included in the Supporting Information, and the time history of the various ice load scenarios, for the location of maximum load change, is shown in Fig. S1. Time histories of the ICE-5 G (Peltier 2004) and ICE-6 G (Argus et al. 2014) models are also shown in Fig. S1 for comparison.

GIA modelling
The ice loading histories described in Section 2.1 are input to a GIA model to solve the sea level equation (Farrell & Clark 1976) and hence determine present-day rates of solid Earth deformation across the Siple Coast. The GIA model accounts for shoreline migration and rotational feedback (Milne & Mitrovica 1998;Mitrovica & Milne 2003;Kendall et al. 2005;Mitrovica et al. 2005), as well as carefully accounting for surface load changes in regions of retreating marine-grounded ice. The latter is achieved by ensuring that retreating marine-grounded ice is immediately replaced by ocean water, where the depth of the water is determined by solving the sea level equation.
There are two components that must be defined in a GIA model; the ice loading history and the rheological properties of the Earth. The ice loading histories described above were interpolated onto a spherical harmonic grid of degree and order 256, corresponding to approximately 80 km spatial resolution. This model resolution is sufficient to investigate the solid Earth response to the modelled ice load changes (predicted uplift rates differ by <0.6 mm yr -1 when the spatial resolution is doubled), and it enables us to combine and compare our results with the previously-published W12 GIA model (Whitehouse et al. 2012b).
The rheological properties of the Earth beneath Antarctica are poorly constrained, but West Antarctica is typically considered to have a thin lithosphere and low viscosity upper mantle compared with East Antarctica (Morelli & Danesi 2004;Lloyd et al. 2013;An et al. 2015). This is supported by a number of studies that have used seismic velocity perturbations to determine lateral and vertical variations in mantle viscosity (and in some cases lithospheric thickness) when deriving a GIA model for Antarctica (Kaufmann et al. 2005;van der Wal et al. 2015). All three studies found upper mantle viscosities to be lower in West Antarctica than East Antarctica, but the magnitude of the difference, and the absolute value of upper mantle viscosity in West Antarctica, differed between studies due to the assumptions made when converting seismic velocity perturbations to mantle viscosities.
Due to the uncertainty in Earth structure, a range of earth models has been used to model the solid Earth response to ice loading. The lithospheric thicknesses used in the modelling are 46, 71 and 96 km and the range of upper mantle viscosities tested is 0.5-5 × 10 20 Pa s. The earth model corresponding to the W12 GIA model (Whitehouse et al. 2012b) is also tested, which has a 120-km-thick lithosphere and an upper mantle viscosity of 1 × 10 21 Pa s. The lower mantle viscosity is fixed at 1 × 10 22 Pa s in all cases.
Although the W12 deglacial history has previously been combined with a specific earth model in order to ensure the best fit between model results and observational data (Whitehouse et al. 2012b), the deglacial model was developed independently of the earth model. This means that it can be combined with other iceloading scenarios, and the response modelled with different Earth parameters, but the resulting model predictions may no longer fit the Antarctic-wide relative sea level data that were used to tune the original W12 GIA model. The results presented here are, therefore, only applicable on a regional scale.

Time stepping
All models include the most recent KIS stagnation event, and for this period the GIA model was run using 10-yr time steps (Table S1). Ice load changes were modelled up to 2010, and for the final 4 yr of each model run 1-yr time steps were used (2010-2014), during which time no load changes occurred. Uplift rates relating to 2012 (the approximate mid-point of the GPS time-series considered in model-data comparisons) are calculated by considering the viscous deformation that takes place between 2010 and 2014; the lack of load changes during this period ensures that there is no elastic component to the model predictions, which could be misinterpreted as GIA.
The KIS-only GIA modelling includes a spin-up period between 3000 and 500 yr BP, during which time it is assumed there were no changes in ice thickness. All other models include load changes prior to 2 ka BP; for the W12+Siple model the time steps from the original W12 model were used (Whitehouse et al. 2012b); for the W12_5k+Siple model the original W12 time steps were again used, but earlier grounding line retreat was achieved by replacing the 5, 4 and 3 ka BP W12 ice sheet configurations with the 2 ka BP W12 ice sheet configuration (i.e. the present-day configuration of the ice sheet). The timings of the post-2 ka BP ice load changes in the W12+Siple and W12_5k+Siple models are given in Table 3.

Accumulation data
Surface mass balance (SMB) data were used to determine the distribution and magnitude of ice thickening along KIS following the completion of stagnation at each location (Section 2.1.1.4). SMB data are available from the regional atmospheric climate model RACMO2.1/ANT (Lenaerts et al. 2012), which provides monthly SMB values between January 1979 and December 2010 on a ∼27 km grid (Fig. 3a). Average SMB rates for this epoch were used when calculating the recent ice loading history of the Siple Coast. Departures from these average rates are evident in local ice core records (e.g. Kaspari et al. 2004;Banta et al. 2008) but they will have had a negligible impact on net ice thickness change over the period of interest. We assumed that the magnitude and spatial pattern of the accumulation rate has been constant for the whole period of interest, and that large fluctuations in accumulation did not occur. These assumptions are justified by noting that RACMO2.1/ANT model output suggests relatively constant accumulation rates over the last ∼30 yr (Fig. 3b). Moreover, centennial accumulation rates deduced from ice cores in the region demonstrate that long-term average accumulation rates are within 0.1 m yr -1 of the RACMO2.1/ANT SMB values (Kaspari et al. 2004;see Fig. 3a).
SMB rates were not used to determine regional ice thickness changes during the earlier episodes of stagnation and reactivation described in Section 2.1.2 since only localised thickness changes were accounted for in these scenarios. However, we note that regional trends in accumulation rates, such as the inferred 0.01-0.05 m yr -1 decrease in accumulation rates at Siple Dome and WAIS Divide over the last 2 ka (e.g. Price et al. 2007;Neumann et al. 2008;Fudge et al. 2013, Supporting Information), could have influenced both the ice dynamics and the ice-loading history of the region.

Ice surface elevation change data
In addition to the SMB data, surface elevation change data were used to constrain reconstructions of recent ice thickness change across the most rapidly thickening portions of the Siple Coast (Section 2.1.1.4). Data are available from ICESat laser altimetry measurements between 2003 February and 2007 November (Pritchard et al. 2009) and CryoSat-2 radar altimetry measurements between 2010 November and 2013 September (McMillan et al. 2014). Together, these data sets show thickening of up to 0.6 m yr -1 on the trunk of KIS during the last decade (Pritchard et al. 2009;Helm et al. 2014). The pattern of thickening is very similar in the two data sets, and the ICESat data were used to constrain our modelling.

GPS data
The results of the GIA modelling were compared with GPSobserved present-day uplift rates at five sites in the region. Details for the GPS sites are given in Table 4, and locations shown on Fig. 1. The observed uplift rates and uncertainties have been taken from Argus et al. (2014). Also given in Table 4 is the modelled elastic uplift at each site based on present-day ice loss in the northern Antarctic Peninsula and Amundsen Sea sector, as taken from Argus et al. (2014). It is worth noting that Gunter et al. (2014) estimate elastic uplift to be <0.3 mm yr -1 at all GPS locations outside of the northern Antarctic Peninsula and Amundsen Sea sector, although they do not include SDLY or RAMG in their analysis. SDLY is located closer to the Amundsen Sea sector than the other sites and is therefore the only site likely to be materially affected by elastic uplift. The nearby site of Whitmore Mountains (WHTM) has been excluded as the time-series is highly non-linear.

Ice loading reconstructions
The results of the ice thickness change calculations relating to the stagnation of KIS are shown in Figs 4-6 and the maximum ice thickness change for each model is given in Table 5. Cumulative change in ice thickness at 10-yr time intervals from 1820 to 2010 is shown for the two end member models, 'LB_100' (Fig. 4) and 'UB_300' (Fig. 6), along with best-estimate model, 'BE_200' (Fig. 5).
While the spatial pattern of thickening is common across all models, the amount of thickening varies. Earlier stagnation generally results in greater total thickening. Similarly, a higher pre-stagnation velocity results in a larger amount of influx and hence increased ice build-up during stagnation.
In all models, a small area of ice build-up, ∼100 km upstream of the main region of thickening, appears in the mid twentieth century. This arises because ICESat elevation change rates are >0.3 m yr -1 in this region (see Fig. 2), and therefore it is included in the ice thickness change calculations once the ice stream has stagnated. This small anomaly will have a negligible effect on the GIA results.
There is significant overlap of the results from the nine ice histories. A shorter stagnation period combined with higher velocity gives similar results to a longer stagnation period with a lower velocity. For example, the maximum ice thickness change from model LB_200 (88.7 m) is similar to that of model BE_100 (81.9 m). Therefore, for input to the GIA model, only three ice histories were used; the extreme upper and lower bounds (LB_100, UB_300), and the best estimate (BE_200). The results for these models are described in the next section. The average maximum ice thickness   change for all models is 100.5 m, almost identical to the 100.4 m from the BE_200 model.

KIS stagnation (KIS-only model)
Model-predicted vertical motion rates due to the build-up of ice following the stagnation of KIS are shown in Fig. 7 for the BE_200 model, for a range of earth models. As expected, the largest subsidence rates occur at the location of maximum ice build-up. For the UB_300 and LB_100 ice models, the amount of subsidence was found to be slightly more or less than the BE_200 model, respectively, for a given earth model. The maximum subsidence for each combination of ice and earth model is given in Table 6. There is potentially a large amount of subsidence related to thickening on KIS, depending on the underlying Earth rheology, with modelled values ranging from -17 mm yr -1 for the weakest earth model combined with largest ice load to around -1 mm yr -1 for the stronger earth models tested. Adoption of the earth model used by Whitehouse et al. (2012b) yields maximum subsidence rates of only -0.2 to -0.3 mm yr -1 (see Table 6 and lower left-hand panel of Fig. 7).

Regional loading model (W12+Siple model)
Predicted present-day uplift rates for the W12+Siple model are shown in Fig. 8. For the weaker earth models, uplift rates are generally in the range -2 to +2 mm yr -1 , apart from across the region of KIS-thickening, where rates reach -7 to -14 mm yr -1 (subsidence). For the stronger earth models the response to Late Holocene load changes is muted and the regional pattern largely resembles the original W12 model (see the last panel of Fig. 8), with a significant amount of uplift predicted beneath the Ross Ice Shelf. For the weaker earth models, particularly those with a thinner lithosphere, uplift is predicted in the region of WIS and MacIS, indicating that the net response to the modelled stagnation-reactivation cycle on these ice streams is present-day rebound. The maximum magnitude of this response is ∼3 mm yr -1 ; much lower than the response to the more recent load changes on KIS. Comparisons between experiments that do and do not include the loading cycle on WIS and MacIS (results not shown) suggest that the ongoing regional uplift associated with ice loss on these ice streams following reactivation acts to damp the subsidence associated with recent KIS thickening. In addition, uplift rates in the region of the present-day grounding line are decreased by up to 1 mm yr -1 when the load-ing cycle on WIS and MacIS is included, due to peripheral bulge collapse following unloading on WIS and MacIS.
The nearest GPS sites are located more than 400 km from the KIS basin, and consequently they cannot be used to determine which model best represents the present-day solid Earth response to KIS stagnation. However, the GPS-derived rates may be used to verify the wider-scale pattern of uplift predicted by the combined W12+Siple model. In general, earth models with an upper mantle viscosity of 1 × 10 20 Pa s (middle column of Fig. 8) show a reasonably good fit to the GPS data, although not appreciably better than the original W12 earth model (last panel of Fig. 8).

Adjusted W12 model (W12_5k+Siple model)
The difference between models that adopt standard and earlier grounding line retreat during the mid-to-late Holocene are shown in Fig. 9 (earlier retreat minus standard W12-modelled retreat). In all cases the earlier-retreat model predicts lower rates of uplift in the region of the Siple Coast grounding line (green areas), and lower rates of subsidence in central West Antarctica (brown areas), with differences peaking at ∼5 mm yr -1 for earth models with an upper  a Models correspond to Fig. 4 (LB_100), Fig. 5 (BE_200) and Fig. 6 (UB_300).
mantle viscosity of 1 × 10 20 Pa s. The difference pattern reflects the more advanced state of relaxation of the solid Earth in the case where post-Last Glacial Maximum (LGM) ice mass changes are not assumed to continue into the Late Holocene. Adopting the earlierretreat model reduces the magnitude of subsidence rates associated with KIS stagnation by several mm yr -1 .

Comparison with empirical GIA solution
For comparison, predicted uplift rates for the two long-term loading histories (W12+Siple and W12_5k+Siple) were compared with the empirically-derived GIA-related uplift rates of Gunter et al.  Table 6. Maximum present-day subsidence rates in response to the KIS-only stagnation model, for the different combinations of ice and earth models. For the earth models, L is lithospheric thickness (km), and UM and LM are upper and lower mantle viscosity (×10 21 Pa s), respectively. Ice history models are described in the main text.
Maximum subsidence (mm yr -1 ) for ice models: In order to compare results from the two approaches, our GIA model results were first smoothed with a 400 km (half width) Gaussian filter, to ensure the resolution of the two solutions was the same, and then differenced with the empirical solution. The difference in uplift rates is shown in Fig. 10, where the two ice loading histories -W12+Siple and W12_5k+Siple-are combined with three different earth models-a weak earth model ( The results for the weaker earth models (first and second columns of Fig. 10) show widespread under-prediction of the Gunter et al. (2014) empirical rates, by up to 10 mm yr -1 . The differences are dominated by subsidence related to the recent build-up of ice on KIS, although a smaller misfit is achieved when post-LGM ice load changes are assumed to be complete by 5 ka BP (second row of Fig. 10). Over the Ross Ice Shelf, there is a better match between the modelled and empirical solutions, with misfits of 3 mm yr -1 or less for the weaker earth models, although we note that this misfit is still greater than the uncertainty on the empirical solution.
For the W12 earth model (third column of Fig. 10), the misfit with the Gunter et al. (2014) model is dominated by large negative values across the Ross Ice Shelf. This indicates that the W12 model and its derivatives over-predict uplift by up to 6 mm yr -1 compared with the Gunter et al. (2014) model when combined with a relatively high upper mantle viscosity. King et al. (2012) suggested that the W12 model over-predicts uplift across the Ross Ice Shelf, estimating that area-averaged rates here should be close to zero. We find the closest agreement with the empirical rates of Gunter et al. (2014) when post-LGM ice load changes are assumed to be complete by 5 ka BP and the local upper mantle viscosity is taken to be 1 × 10 20 Pa s (lower middle plot, Fig. 10). However, based on comparison with the empirical solution, we note that the response to KIS thickening seems to be over-predicted when this earth model is adopted.
Differences between GPS-observed uplift rates and GIAmodelled uplift rates (this study) are also shown on Fig. 10. GPS misfit values are generally different to the empirical solution misfit values, which is to be expected since the empirical solution does not agree well with the GPS uplift rates in this region (Gunter et al. 2014). Misfits between GPS uplift rates and GIA model predictions are ≤±1 mm yr -1 at SDLY and HOWE for all combinations of ice and earth models. Misfits at RAMG are up to ±3 mm yr -1 , although this discrepancy is smaller than the observational uncertainty (±5.1 mm yr -1 ; see Table 4); this is the only site to agree with the empirical solution. The largest misfits are found at CLRK and PATN, where the GIA models under-predict uplift by up to 4.5 mm yr -1 . Estimated elastic uplift due to present-day ice loss at these sites is only ∼0.4-0.6 mm yr -1 (Table 4), which cannot explain the misfit, however, we do note again that observational uncertainties are relatively large at these sites (Table 4).

Ice loading history reconstruction
Using the available observational data a set of plausible ice loading histories have been constructed for the ice build-up relating to the stagnation of KIS. Uncertainties in the timing of ice stream stagnation, as given by Retzlaff & Bentley (1993), have been explored to obtain lower and upper bound ice histories. The pre-stagnation velocity of KIS, which affects how much influx and build-up occurs, is also unknown and three realistic velocities were tested. This results in a total of nine ice loading histories which produce realistic amounts of ice thickness change, between 70 and 130 m over a time span of 148-188 yr (see Table 5).
In determining these bounds on net ice thickening we adopted the timings detailed in Retzlaff & Bentley (1993). However, Shabtaie & Bentley (1987) estimate the burial of crevasses, and hence the timing of stagnation, to be around 100 yr before that estimated by Retzlaff & Bentley (1993). Earlier stagnation would result in greater ice    thickness increase and potentially greater present-day subsidence. If this were the case, the results presented in this paper represent a lower bound for the effect of recent ice-load change on solid Earth deformation across the Siple Coast.
One of the limitations of the method used to calculate ice build-up is the assumption that stagnation of a given grid cell occurs in 1 yr. Stagnation at a given location likely occurred over several years, as indicated by the Duckfoot region that shutdown in around 10 yr (Catania et al. 2006). The result of this assumption is that ice thickness change would be underestimated: if the ice stream shut down over a number of years there would be additional influx to the system which has not been accounted for in these calculations. However, the amount of underestimation is likely to be small and within the bounds already covered by varying the timing and velocity, which have the biggest impact on ice thickness change. Furthermore, the effect of this on the resulting GIA signal would be negligible.
The ice loading histories have been constructed using a simple approach which does not incorporate modelling of ice-sheet dynamics. In this region, dynamic ice-sheet models struggle to accurately reproduce ice streams, and models have been specifically adapted to mimic observed behaviour by altering basal sliding conditions (Whitehouse et al. 2012a). Ice-sheet models using a shallow ice approximation are not capable of fully reproducing the behaviour of ice streams (Kirchner et al. 2011). Employing a more sophisticated modelling technique, Pollard & DeConto (2012) derived basal sliding iteratively by matching modelled ice-sheet surface elevations to observed ice-sheet surface elevations. However, this method failed to reproduce observed velocities for KIS. In the future, a more detailed reconstruction of the Siple Coast is needed, but will require a more sophisticated ice-sheet model to successfully capture the ice stream behaviour (e.g. Golledge et al. 2014). For the purposes of this study, our simplified approach is sufficient to enable the sensitivity of GIA to Late Holocene perturbations to be investigated.

Ice loading history data
Ice thicknesses are needed to calculate ice influx (eq. 1), and values were taken from Bedmap2 (Fretwell et al. 2013), which has been compiled from data collected over the past ∼50 yr. However, using near-present-day ice thickness values may lead to an overestimation of net ice-thickness increase because the ice sheet would have been thinner prior to stagnation and hence ice flux into the stagnated area would have been less. The effect of this was tested by subtracting the net change in ice thickness for each model from the initial Bedmap2 ice thicknesses, yielding an estimate for the ice-sheet thickness prior to the stagnation of KIS. The ice histories were calculated again with the adjusted Bedmap2 ice thickness and it was found that the maximum ice-thickness change was up to 5 m less. This is equivalent to a <5 per cent difference compared with the original calculations, which has a negligible effect on the GIA model estimates.
Accumulation rates play a role in reconstructing the ice loading history on KIS because precipitation contributes to net ice thickening downstream of stagnation, where there is no flow to balance mass input to the system. We used accumulation rates from RACMO2.1/ANT, and it was assumed that the magnitude and spatial pattern of accumulation remained constant for the whole time period of the stagnation. Any errors due to this assumption would be small and this was tested by calculating the ice loading histories again with double and half the accumulation rate. Differences to the maximum ice-thickness change for each model were <5 per cent, which is well within the bounds captured by our end-member models.

GIA due to KIS stagnation
When modelling the build-up of ice related to the stagnation of KIS, the GIA model-predicted subsidence varies depending on the combination of ice and earth model used (see Table 6). For the best estimate ice model with a maximum ice thickness increase of 100 m, subsidence is in the range -9 to -14 mm yr -1 for an upper mantle viscosity of 5 × 10 19 Pa s, which reduces to -5 to -7 mm yr -1 for 1 × 10 20 Pa s. For an upper mantle viscosity of 5 × 10 20 Pa s and higher (e.g. the W12 earth model), subsidence rates are negligible. The magnitude of subsidence increases for the ice loading histories with larger amounts of ice-build-up (UB_300 ice model) and decreases with lower amounts of ice build-up (LB_100 ice model). This demonstrates that there may be a significant present-day GIA signal from stagnation-related ice build-up over the past 165 yr, although it heavily depends on the regional upper mantle viscosity. Our results suggest that this recent ice history should be included in ice loading models if the upper mantle in this region is less than 5 × 10 20 Pa s.

GIA on the Siple Coast
Ice build-up due to the recent stagnation of KIS has been considered in the context of a post-LGM deglacial model by combining the KIS ice loading history with a model of regional ice load change across Siple Coast during the last 2000 yr, and the W12 deglacial model. A number of different rheological models were considered, and a caveat to this approach is that only local GIA effects can be considered when the W12 deglacial history is combined with earth models that are different to the W12 best-fitting earth model.
For weaker earth models, much of the uplift signal due to largescale changes since the LGM has diminished by the present-day, and the dominant signal is due to loading on KIS (Fig. 8). For an upper mantle viscosity of 5 × 10 20 Pa s and higher, the magnitude and spatial pattern of uplift is similar to that of the original W12 uplift as the recent loading on the Siple Coast has very little effect (<1 mm yr -1 ) on the present-day uplift. Comparison of GIA modelpredicted uplift with GPS-observed uplift (Fig. 8) provides limited insight because the GPS sites are not located close enough to KIS to constrain deformation due to stagnation-related ice build-up.
It is important to include all post-2 ka BP ice load changes across the wider Siple Coast because changes over these time scales can play a dominant role in the present-day rebound signal if the regional upper mantle viscosity is weak. Evidence exists for multiple stagnation and reactivation cycles on the Siple Coast ice streams (Hulbe & Fahnestock 2007;Catania et al. 2012), which would cause fluctuations in ice thickness as ice builds up during stagnation and thins on reactivation. We model just a single cycle of loading and unloading on each of KIS, WIS and MacIS prior to the most recent stagnation of KIS, and our results suggest that for a weak earth model these earlier load changes play a role in defining the pattern and magnitude of present-day rebound: uplift due to ice-mass loss associated with the earlier reactivation of the Siple Coast ice streams acts to damp the significant subsidence from the most recent loading on KIS, while rebound rates in the southern Ross Sea are damped by up to 1 mm yr -1 due to peripheral bulge subsidence to the north of the reactivated ice streams. Further investigation into the detailed history of ice load changes on the Siple Coast is warranted as it has been demonstrated that the response to cyclic loading depends on the magnitude and timing of the loading cycles (Ivins et al. 2000).

Dependence on Holocene ice load model
Between 5 and 2 ka BP the W12 deglacial model prescribes grounding line retreat and ice thinning in the region of the present-day Siple Coast grounding line, and a small amount of ice thickening across central West Antarctica in response to warming temperatures and hence increasing accumulation rates (Muszynski & Birchfield 1985;Siegert & Payne 2004). The response to these load changes dominates the uplift pattern when the original W12 model is combined with a weak earth model (Whitehouse et al. 2012b), leading to a spatial pattern of deformation that is very different to that predicted by Gunter et al. (2014). Assuming the results of Gunter et al. (2014) are robust, this mismatch may be due to (i) the neglect of Late Holocene load changes along the Siple Coast, (ii) the use of an unfeasibly weak earth model or (iii) the incorrect modelling of mid-Holocene ice load changes in this region. The first point is the main focus of this study, while the second is an important issue, which is hard to address due to the lack of independent constraints on mantle rheology and the difficulty of directly measuring solid Earth rebound in this region. The third point is explored in our final set of experiments (Section 3.4) and discussed further below.
Radio-echo sounding data suggest that the majority of the transition from glacial to interglacial accumulation rates was complete by ∼6.4 ka BP (Siegert & Payne 2004), hence ice thickening may have been minimal during the Late Holocene. There is also evidence to suggest that grounding line retreat in the central Ross Sea may have been faster than the retreat recorded along the western margin of the Ross Sea (Anderson et al. 2014 and references therein;McKay et al. 2015), and hence largely complete by 5 ka BP. Both of these issues were addressed here by altering the W12 deglacial model such that large-scale ice load changes are complete by 5 ka BP, and the only load changes subsequent to this occur during the last 2000 yr along the Siple Coast. The W12 model is the product of a time slice approach rather than a transient ice sheet reconstruction, and therefore there is flexibility in how ice thickness changes are interpolated between the main temporal reconstructions, that is between 5 ka BP and present (Whitehouse et al. 2012a).
The use of a model which assumes that large scale, post-LGM ice sheet changes are complete by 5 ka BP, and that subsequent changes only occur along the Siple Coast, results in lower magnitude predictions for present-day uplift and subsidence rates, and an improved agreement with empirically-derived models of at University of Durham on March 31, 2016 http://gji.oxfordjournals.org/ Downloaded from GIA (Fig. 10). However, the non-uniqueness of the problem implies that further work is required to independently constrain both the Earth rheology of this region and the ice sheet history.

Comparison with observations
Results from our forward GIA modelling, for three earth models (weak, medium and the W12 earth model; see Fig. 10), were compared with the Gunter et al. (2014) empirical GIA estimate. This estimate is derived from GRACE and ICESat observations and is independent of ice or earth model assumptions, but is sensitive to data analysis approaches and the accuracy of a firn densification model.
Differences between the empirical estimate and GIA models that adopt a relatively weak earth model are dominated by relative subsidence at KIS (indicated by the red areas on Fig. 10, left-hand and centre columns). The absence of this subsidence signal in the Gunter et al. (2014) estimate may indicate that the mantle viscosity in this region is relatively high, greater than 1 × 10 20 Pa s. The misfit is significantly reduced, by up to 4 mm yr -1 , when post-LGM ice load changes are assumed to be complete by 5 ka BP (bottom row of Fig. 10).
When the W12 earth model is adopted, the misfit with the empirical estimate is dominated by the large uplift centre over the Ross Ice Shelf that is present in the original W12 model (Fig. 10, right-hand column). Differences peak at 6 mm yr -1 ; this is larger than the 1σ uncertainty on the empirical solution, which peaks at ∼2 mm yr -1 in West Antarctica (Gunter et al. 2014). The misfit over the Ross Ice Shelf is reduced to <3 mm yr -1 for the weaker earth models (Fig. 10, centre and left-hand panels); however, it is not necessarily the case that this improved agreement implies the Earth must be weak because the misfit may also be due to errors in the ice history. We have explored uncertainties associated with the timing of ice load changes in the W12 deglacial model, but if, for example, the W12 model assumes too much ice loss between the LGM and present, this could also explain the over-prediction of present-day uplift rates. Consequently, the cause of the misfit is not clear. Comparison with the empirical GIA solution of Schoen et al. (2015) would reveal larger misfits as they predict consistently lower uplift rates than Gunter et al. (2014). Misfits between GIA model-predicted uplift and the GPSobserved uplift, shown in Fig. 10, are within ±1 mm yr -1 for two of the five sites considered, and within ±3 mm yr -1 for a third site. The largest misfits are found at CLRK and PATN, where our models under-predict uplift by up to 4.5 mm yr -1 . However, the uncertainties on uplift rates at these sites are large (Table 4) and the differences are not significant. The nearest GPS sites (Argus et al. 2014) are located far from the main areas of misfit, and uplift rates at these sites have large uncertainties, meaning that they cannot support or disprove the main conclusions drawn from comparison with the empirical solution. The absence of rock outcrops in the central Siple Coast (Fig. 10) means that GPS uplift rates cannot play a substantial role in separating competing GIA model predictions for this region.

C O N C L U S I O N S
(1) A simple ice loading history has been constructed for the recent stagnation of KIS resulting in between 70 and 130 m ice thickness increase over a period of 148-188 yr.
(2) GIA model-predicted subsidence in response to the stagnation of KIS peaks at -17 mm yr -1 for the weakest earth models but is much less for models with a higher mantle viscosity (-4 to -8 mm yr -1 for 1 × 10 20 Pa s and less than -1 mm yr -1 for a viscosity of 5 × 10 20 Pa s).
(3) Combining ice load fluctuations over the last 2000 yr along the Siple Coast with the W12 deglacial model shows that, for a weak earth model, much of the LGM signal has diminished and present-day uplift is dominated by subsidence over KIS, although KIS-related subsidence is partly cancelled out by uplift associated with ice unloading several hundred years earlier. For the stronger earth models, the post-2 ka BP load changes explored in this study have a negligible impact on present-day uplift rates.
(4) Comparing the GIA model predictions with an empirical GIA estimate (Gunter et al. 2014) produces misfits that are dominated by subsidence across KIS for the weaker earth models, and uplift across the Ross Ice Shelf for the stronger earth models. The absence of any large subsidence in the empirical estimate suggests that the upper mantle viscosity beneath the Siple Coast is greater than 1 × 10 20 Pa s. Misfits are likely due to a combination of errors in the ice history and earth model.
(5) Adopting an earlier retreat for the W12_5k+Siple model alters the magnitude of subsidence due to the stagnation of KIS by several mm yr -1 , and when combined with an upper mantle viscosity of 1 × 10 20 Pa s it improves the fit with the Gunter et al. (2014) empirical GIA solution.
(6) Our revised predictions of GIA uplift rates are small over the Ross Ice Shelf and Siple Coast only when considering upper mantle viscosities of 0.5-1.0 × 10 20 Pa s (Fig. 8). Such models would be compatible with the finding of King et al. (2012)-that area-averaged Ross Ice Shelf GIA is near zero-whilst still allowing substantial ice mass loss since the LGM.
(7) Late Holocene ice load changes related to the stagnation and reactivation of ice streams on the Siple Coast may play a dominant role in defining the pattern of present-day uplift. In the future, numerical ice-sheet models that are able to capture dynamic ice stream behaviour should be used to determine a more detailed loading history for this region.

A C K N O W L E D G E M E N T S
GAN was supported by a NERC PhD Studentship, PLW is a recipient of a NERC Independent Research Fellowship (grant number NE/K009958/1) and MAK is supported by Australian Research Council Future Fellowship (project number FT110100207) and Special Research Initiative for Antarctic Gateway Partnership (Project ID SR140300001). We are grateful to Hamish Pritchard for providing the ICESat data shown in Fig. 2 and Brian Gunter for providing model output used to produce Fig. 10. We are also grateful to Donald Argus and two anonymous reviewers for providing constructive reviews of the manuscript.