The impact of lateral variations in lithospheric thickness on glacial isostatic adjustment in West Antarctica

Differences in predictions of Glacial Isostatic Adjustment (GIA) for Antarctica persist due to uncertainties in deglacial history and Earth rheology. The Earth models adopted in many GIA studies are deﬁned by parameters that vary in the radial direction only and represent a global average Earth structure (referred to as 1-D Earth models). Oversimplifying the actual Earth structure leads to bias in model predictions in regions where Earth parameters differ signiﬁcantly from the global average, such as West Antarctica. We investigate the impact of lateral variations in lithospheric thickness on GIA in Antarctica by carrying out two experiments that use different rheological approaches to deﬁne 3-D Earth models that include spatial variations in lithospheric thickness. The ﬁrst experiment deﬁnes an elastic lithosphere with spatial variations in thickness inferred from seismic studies. We compare the results from this 3-D model with results derived from a 1-D Earth model that has a uniform lithospheric thickness deﬁned as the average of the 3-D lithospheric thickness. Irrespective of the deglacial history and sublithospheric mantle viscosity, we ﬁnd higher gradients of present-day uplift rates (i.e. higher amplitude and shorter wavelength) in West Antarctica when using the 3-D models, due to the thinner-than-1-D-average lithosphere prevalent in this region. The second experiment uses seismically inferred temperature as an input to a power-law rheology, thereby allowing the lithosphere to have a viscosity structure. Modelling the lithosphere with a power-law rheology results in a behaviour that is equivalent to a thinner lithosphere model, and it leads to higher amplitude and shorter wavelength deformation compared with the ﬁrst experiment. We conclude that neglecting spatial variations in lithospheric thickness in GIA models will result in predictions of peak uplift and subsidence that are biased low in West Antarctica. This has important implications for ice-sheet modelling studies as the steeper gradients of uplift predicted from the more realistic 3-D model may promote stability in marine-grounded regions of West Antarctica. Including lateral variations in lithospheric thickness, at least to the level of considering West and East Antarctica separately, is important for capturing short-wavelength deformation and it has the potential to provide a better ﬁt to Global Positioning System observations as well as an improved GIA correction for the Gravity Recovery and Climate Experiment data.

component of the total gravitational signal and must be removed to yield estimates for ice-mass balance .
Traditionally, global and Antarctic-wide models of GIA have used a 1-D approximation of the Earth structure consisting of an elastic lithosphere underlain by a linear viscoelastic upper and lower mantle, where properties vary only radially (e.g. Peltier 1974;Milne & Mitrovica 1996;Kendall et al. 2005). In reality, the structure of the Earth is far more complex and models that reflect lateral as well as vertical variations in Earth properties are needed to provide more accurate predictions of present-day GIA-related deformation and geoid changes, both in Antarctica van der Wal et al. 2015;Sasgen et al. 2017) and elsewhere, for example, Greenland (Khan et al. 2016). Including 3-D structure in GIA models is particularly pertinent for Antarctica as this continent is considered to consist of two distinct regions in terms of Earth structure: a thick cratonic lithosphere and a high-viscosity uppermost mantle in the East, and thinner lithosphere and lower viscosity uppermost mantle in the West (Morelli & Danesi 2004). Modelling East and West Antarctica with a 1-D Earth model as described above, therefore, has the potential to produce incorrect estimates of the presentday GIA signal in one or both of these regions. For example, A et al. (2013) compared deformation rates predicted by a 3-D model incorporating laterally varying lithospheric thickness and mantle viscosity with a model that is the 1-D average of the 3-D profile and found mismatches at Global Positioning System (GPS) locations in Antarctica. Furthermore, capturing variability in the Earth structure within West Antarctica is important because regional 1-D GIA studies have indicated differences in the Earth structure across the region (Nield et al. 2014;Wolstencroft et al. 2015;Nield et al. 2016;Zhao et al. 2017).
This study focusses on how lateral variations in lithospheric thickness impact predictions made by GIA models. The lithospheric thickness can be defined by various criteria, such as a change in the method of heat transfer (Martinec & Wolf 2005), seismic anisotropy or resistivity (Eaton et al. 2009). In GIA modelling, the lithosphere is defined on the basis of mechanical properties and is considered to be a part of the crust and upper mantle that behaves elastically on timescales of glacial cycles (tens of thousands of years; Martinec & Wolf 2005;Watts et al. 2013;Kuchar & Milne 2015). The lithosphere can be modelled with either a purely elastic rheology, that is, it has no viscous component (e.g. Argus et al. 2014), or as a viscoelastic material with sufficiently high viscosity that it does not relax in response to surface loading on timescales of a glacial cycle (e.g. Kendall et al. 2005;Whitehouse et al. 2012b;Kuchar & Milne 2015), thereby behaving elastically. Studies have also combined these approaches, for example, Kaufmann et al. (2005) modelled a 100 km thick lithosphere composed of a 30 km purely elastic layer overlying a 70 km viscoelastic layer with a viscosity of 1 × 10 24 Pa s, which is approximately the limit of what could be considered elastic over GIA timescales [e.g. 1 × 10 22 Pa s (Sasgen et al. 2017) -1 × 10 24 Pa s (Khan et al. 2016)]. Kuchar & Milne (2015) investigated the effect of depth-dependent viscosity in the lithosphere on relative sea-level predictions using a radially varying (i.e. 1-D) Earth model and found that predictions made using a lithosphere with viscosity structure were similar to predictions made using a purely elastic, but much thinner, lithosphere.
To some extent, the apparent thickness of the lithosphere depends on the timescales of surface loading. Over long timescales (∼1 Myr), viscous relaxation in the lower lithosphere means that the lithosphere seems to behave as a relatively thin elastic layer (Watts et al. 2013). However, over GIA timescales (∼100 kyr), the lithosphere seems to behave as a thicker elastic layer (Martinec & Wolf 2005;Nield et al. 2014;Wolstencroft et al. 2014). On the basis of wave speed variations, seismic studies can distinguish between thermal conduction and convection regimes in the upper mantle. The conductive domain defines the tectonic plate. However, the elastic thickness varies as a function of timescale of surface loading and is typically thinner than the seismic lithospheric thickness.
The studies and methods described above have used linear viscoelastic rheology to model GIA. However, the use of power-law rheology is becoming increasingly common (Wu 1999;Barnhoorn et al. 2011;van der Wal et al. 2013;van der Wal et al. 2015). van der Wal et al. (2015) used seismic velocity anomalies (Grand 2002) and geothermal heat flux estimates (Shapiro & Ritzwoller 2004) for Antarctica to infer mantle temperatures that were used to derive creep parameters for input to a power-law rheology. By defining spatially varying creep parameters, the GIA model included laterally varying Earth structure. For this approach, the lithospheric thickness is implicitly defined by the creep parameters, rheological model and some threshold viscosity above which it can be considered to behave elastically as described above.
Modelling advances in the past few decades (Wu & Johnston 1998;Latychev et al. 2005b;van der Wal et al. 2013) have eased the computational burden of 3-D GIA modelling and detailed data sets are now available that can be used to define lateral Earth structure (Ritsema et al. 2011;Heeszel et al. 2016). Hence, there are an increased number of studies incorporating 3-D Earth structure into GIA models with both linear and nonlinear rheologies. Several approaches can be used to infer 3-D mantle viscosity (Ivins & Sammis 1995;Kaufmann et al. 2005) and lithospheric thickness for input to GIA models, with the latter being the focus of this study. A seismically derived lithosphere-asthenosphere boundary (LAB) depth is sometimes used to infer laterally varying GIA lithospheric thickness with linear viscoelastic rheology, after scaling to account for differences between a seismically derived definition of the lithosphere and the mechanical definition used in GIA studies. For example, Kaufmann et al. (2000) reduced a seismically derived LAB depth by a factor of two for their GIA modelling, and Khan et al. (2016) used an adjustment factor to scale the LAB depth published by Priestley & McKenzie (2013). However, it is not clear whether a lithosphere defined by seismic properties can be converted to a lithosphere defined by mechanical properties through a scaling factor. Seismic properties could have a different dependence on temperature and composition than mechanical properties. One way to circumvent this issue is to use temperatures derived from seismic velocity perturbations as input to a power-law rheology, which eliminates the need to explicitly define lithospheric thickness (van der Wal et al. 2013). In this approach, assumptions are made in converting seismic velocity anomalies into temperature and viscosity, and the lithosphere is defined implicitly by the effective viscosity. So although temperatures ultimately come from the same seismic source as the LAB depths, no new assumptions are required other than those made for converting seismic velocities to viscosity.
Previous studies investigating 3-D Earth structure in GIA models of Antarctica have focussed on the effect of lateral variations in mantle viscosity (e.g. Kaufmann et al. 2005) or a combination of laterally varying lithospheric thickness and upper-mantle viscosity van der Wal et al. 2015) on present-day uplift rates. Studies that isolate the effect of including lateral variations in lithospheric thickness in models of GIA exist for regions in the northern hemisphere (Kaufmann et al. 2000;Zhong et al. 2003;Latychev et al. 2005a;Whitehouse et al. 2006;Steffen et al. 2014) but currently no such study exists for Antarctica. The aim of this study is to isolate the effect of lateral variations in lithospheric thickness on GIA in West Antarctica to determine the effect on gradients of present-day uplift rates when compared with a 1-D Earth model. We explore the two methods of defining a laterally varying lithospheric thickness mentioned above. The first method (experiment 1) uses a scaled seismically inferred LAB depth to determine spatial variability of an elastic lithosphere. For this method, we employ two different models of seismically derived LAB depth (experiments 1a and 1b). We present results that focus on the differences in gradients of present-day uplift rate between 1-D and 3-D models. The difference in the spatial gradient of uplift rate indicates how the amplitude and wavelength of deformation varies between the two models. Each 3-D model includes lateral variations in lithospheric thickness derived from one of the two seismic models and the equivalent 1-D model has a uniform lithospheric thickness that is simply the average of the lithospheric thickness in the 3-D model. Using this method, we seek to determine to what degree the differences between 1-D and 3-D models are independent of the details of: (1) the assumed deglacial history, and (2) the sublithospheric upper-mantle viscosity.
The second method (experiment 2) uses seismically inferred temperature as input to a power-law relationship, thereby assigning a viscosity structure to the lithosphere. In this method, the lithospheric thickness is implicitly defined as the depth at which the resulting viscosity is small enough so that significant deformation takes place during a glacial cycle. In order to determine the effect of including viscosity in the lithosphere through a power-law rheology on gradients of present-day uplift rate, we compare results using the power-law rheology to those from the first method which assumes that the lithosphere is elastic. For both methods, we use a finite element model with 20 km thick layers, meaning that the variable lithospheric thickness is captured in 20 km 'steps' between element locations. Due to large uncertainties in both Earth structure and ice history, we do not attempt to fit any observational data such as GPS-observed uplift.

Model and geometry
We use a 3-D flat-earth finite element model constructed with the ABAQUS software package (Hibbitt et al. 2016) to compute the solid Earth deformation in response to a changing surface load. The validity of using this approach to model the Earth's response to changes in ice-sheet loading has been shown previously by Wu & Johnston (1998). This method has been used in many studies to model GIA in regions such as Fennoscandia (Kaufmann et al. 2000;Steffen et al. 2006), Antarctica (Kaufmann et al. 2005) and Iceland (Auriac et al. 2013), and has the advantage of computational efficiency over a spherical global model. The flat-earth finite element approach has been shown to be accurate when computing deformation within the ice margin for ice loads with comparable size to the Laurentide Ice Sheet (Wu & Johnston 1998), which makes it applicable to the Antarctic Ice Sheet with its smaller lateral extent.
The mesh consists of eight-node brick elements with reduced integration. The surface geometry of the mesh consists of a 3500 × 4500 km area of interest embedded in a larger model domain. The area of interest represents West Antarctica and has elements of 100 × 100 km (elements are shown in Fig. 1). Outside this region, the element size increases towards the periphery of the model, from 550 km in East Antarctica to approximately 5000 km at the edge of the model domain, for computational efficiency, and the domain has an overall width of 60 000 km. The extremely large model domain is required to ensure that boundary effects are negligible in the area of interest (Steffen et al. 2006). We do not model any ice loading changes outside Antarctica as the impact on uplift rates in Antarctica would be negligible (Whitehouse et al. 2012b) but we do include ocean loading. The model consists of 22 depth layers representing the Earth's surface to the core-mantle boundary (Table 1). A 30 km purely elastic crustal layer [the same thickness used by Kaufmann et al. (2005)] is underlain by eleven 20 km thick layers to 250 km depth to capture a higher resolution in the lithosphere and uppermost mantle. Below this, layers are thicker (i.e. lower resolution with depth) since surface deformation will be less sensitive to the details of mantle rheology below 250 km depth (Lau et al. 2016). The buoyancy force is accounted for by applying Winkler foundations to layer boundaries where a density contrast occurs (Wu 2004). To ensure a direct comparability between the models, the same mesh is used for all the experiments.

Earth models and data
The compressible elastic material properties for each layer described above are listed in Table 1. The elastic and density structures of the Earth are derived from the preliminary reference Earth model (Dziewonski & Anderson 1981). For each element below the uppermost elastic layer, creep parameters are assigned on an element-byelement basis. The geometry of the mesh means that the lithospheric thickness (experiment 1) or the viscosity (experiment 2) can vary in 20 km steps between the adjacent element locations. Latychev et al.  (2005a) demonstrated that the differences in uplift rate over previously glaciated regions in the northern hemisphere peak at 2 mm yr −1 for a jump of 150 km between continental and oceanic lithospheric thickness, so we conclude that the effect of a 20 km jump on the uplift rate is likely to be small. Our approach to defining variable lithospheric thickness allows us to use the same mesh for all experiments, thus ensuring that the results are directly comparable. The sublithospheric upper mantle (down to a depth of 670 km) is a linear viscoelastic layer with uniform viscosity and several different upper-mantle viscosities are tested to determine dependence of the results on the underlying viscosity (see Table2). Properties for the lower mantle are the same for all models (Table 1). The thickness of the lithosphere is defined differently in experiments 1 and 2, as detailed in Table 2. The first experiment considers an elastic lithosphere with spatial variability defined by two different models of seismically derived LAB depth (Priestley & McKenzie 2013;An et al. 2015a), as described in the following sections. Seismically derived LAB depths tend to be thicker than those inferred from GIA studies. For example, Steffen et al. (2014) compare GIA-inferred lithospheric thicknesses in the Baltic Sea with three LAB depth models and find a consistently thinner lithosphere by 30-80 km. We, therefore, uniformly reduce the LAB depths from the Priestley & McKenzie (2013) and An et al. (2015a) models so that the thicknesses are more representative of a GIA lithospheric thickness. We use two models to test if the resolution and accuracy of the seismically derived LAB is important. The An et al. (2015a) model gives a greater level of detail, and hence more spatial variability in lithospheric thickness, than that of Priestley & McKenzie (2013) because it is an Antarctic-specific model derived using many additional seismic stations. For those elements representing the lithosphere, the viscosity is set to 1 × 10 49 Pa s to mimic elastic behaviour on glacial timescales, apart from the uppermost 30 km layer which is modelled as purely elastic. This approach of modelling an elastic lithosphere using a viscoelastic rheology with high viscosity (including combinations of purely elastic and viscoelastic rheology) is taken in many GIA studies (Peltier 2004;Kaufmann et al. 2005;Kendall et al. 2005;Whitehouse et al. 2012b;Ivins et al. 2013;Wolstencroft et al. 2015) and for the timescales we are interested in the lithosphere is generally considered to be elastic at viscosities above 1 × 10 24 -1 × 10 25 Pa s (Kaufmann et al. 2005;Steffen et al. 2006;Barnhoorn et al. 2011;Khan et al. 2016). Throughout the rest of this paper, we simply refer to this type of modelled lithosphere as the 'elastic lithosphere'. For each model in experiment 1, we compare the laterally varying model with an equivalent 1-D model in which the lithospheric thickness is simply the average of the 3-D lithospheric thickness.
The second experiment uses a power-law rheology; this complementary approach allows us to investigate the differences in the two methods used to define the lithospheric thickness in GIA modelling. Mantle temperatures (An et al. 2015a) are used to determine diffusion and dislocation creep parameters following the methods described by Hirth & Kohlstedt (2003) and van der Wal et al. (2013;2015). The reader is referred to these papers for a detailed description of the method. We limit the power-law rheology to the same horizontal and vertical domain as defined by the An et al. (2015a) LAB depths for two reasons: (1) so that the results of experiment 2 can be directly compared with the results of experiment 1, and (2) so that the upper-mantle viscosity remains laterally uniform and, therefore, the effect of a spatially variable lithospheric thickness is isolated from all other parameters.

Priestley & McKenzie (2013)
Priestley & McKenzie (2013) published a model of global seismic velocities from the surface wave tomography which they used to derive mantle temperatures and lithospheric thickness on a 2 • grid with a resolution of 250 km horizontally and 50 km vertically. The lithospheric thickness is defined by the change in heat transfer from conduction to convection. We use this information in experiment 1a to define an elastic lithospheric thickness. The authors state that the uncertainty on the lithospheric thickness is 20-30 km, and therefore   Fig. 1(a) shows the adjusted lithospheric thicknesses mapped onto the ABAQUS layers (Table 1), that is, at each location on the mesh, the adjusted LAB depth is rounded to the nearest layer boundary. Fig. 1(b) shows where the lithosphere in the 3-D model is thicker or thinner than the mean of the LAB depths (90 km, calculated over the region south of 60 • S). West Antarctica has a thinner than average lithosphere whereas East Antarctica has a thicker than average lithosphere.

An et al. (2015a,b)
The second model used in this study is from An et al. (2015a), who infer temperatures below Antarctica from the 3-D seismic velocity model AN1-S (An et al. 2015b), which has a horizontal resolution that increases from ∼120 km in the crust to ∼500 km at a depth of 120 km and a vertical resolution of ∼25 to 150 km depth followed by ∼50 to 250 km depth. We use this model in two ways. First, the seismically derived LAB, which is defined by the depth where the adiabat crosses the 1330 • C geotherm, is used to define lithospheric thickness for the elastic lithosphere in experiment 1b. The uncertainty on the temperature is reported to be ±150 • C, equivalent to ±15-30 km for the LAB depth, so we reduce the LAB depth by 15 km to be representative of GIA-elastic thicknesses, as explained in Section 1, and to be consistent with the scaling of the Priestley & McKenzie (2013) model. Second, we use the temperatures directly as input to the power-law rheology in experiment 2. In this experiment, we consider the 3-D spatial domain that defines the lithosphere in experiment 1, but within this domain we use a power law instead of elastic rheology. Results show the comparison of the two 3-D models, thereby highlighting the differences due to rheological definitions (see Table 2 for a summary of the models). In their model An et al. (2015a) do not infer temperatures for depths shallower than 55 km. Therefore, when using the temperatures in our model we specify a second elastic layer between 30 and 50 km depth to bridge the gap between our uppermost elastic layer and the temperature inputs.
The LAB depths mapped onto the model elements is shown in Fig. 1(c), again with Fig. 1(d) showing where the 3-D model has thicker or thinner lithosphere than the 1-D average [90 km, calculated over the Antarctic Plate which is the spatial limit of the inferred LAB depths in the An et al. (2015a) model]. Similar to Priestley & McKenzie (2013), the lithosphere under East Antarctica is thicker than the 1-D average for the An et al. (2015a) model. The location of the boundary between East and West Antarctica is similar in both seismic models, indicating that the uncertainty on the location of the boundary is small. Some isolated regions of anomalously thick lithosphere are also present in the Northern Antarctic Peninsula, which the authors attribute to a remnant subducted slab from the former subduction zone in this region.

Ice loading
The deglacial history of Antarctica since the LGM is poorly known due to a lack of constraining data and consequently, there remain large differences between recent deglacial models (Whitehouse et al. 2012a;Briggs & Tarasov 2013;Gomez et al. 2013;Ivins et al. 2013;Argus et al. 2014). Given this uncertainty, one of the aims of this study is to investigate whether differences between 1-D and 3-D models are independent of the assumed ice history. To test this, we use several different deglacial models, or ice loading histories, that have quite different spatial patterns and magnitude of loading changes, which, when applied to a specific Earth model, give different patterns of deformation. We compare gradients of present-day uplift rates between 1-D and 3-D Earth models using the same ice history, revealing differences that may be directly attributed to the introduction of a varying lithospheric thickness, and then qualitatively compare the results from different ice models.
Three ice loading scenarios are used in the modelling: W12 (Whitehouse et al. 2012a), ICE-5G (Peltier 2004) and its successor ICE-6G C (Argus et al. 2014;Peltier et al. 2015). The W12 ice loading model was derived using a glaciologically consistent numerical ice-sheet model that was tuned to provide the best possible fit to constraints of ice thickness change, whereas ICE-5G and ICE-6G C have been tuned to fit ice thickness change constraints and GPS-observed uplift rates without satisfying ice-sheet physics. Furthermore, in an attempt to fully isolate differences associated with the introduction of a laterally varying lithospheric thickness from those caused by spatial variations in ice loading, we also construct an idealistic, spatially uniform loading history. In this scenario, the amount of ice thickness change since the LGM is spatially uniform over the grounded area of the present-day ice sheet (as shown in Table3 and Fig. 2)

Ocean loading
The model approach we have used in this study does not solve the sea-level equation (Farrell & Clark 1976) and cannot compute variable sea level with time. We, therefore, take the approach of applying an ocean load that has been derived using a global, spherically symmetric GIA model (Mitrovica & Milne 2003;Kendall et al. 2005;Mitrovica et al. 2005). The GIA model uses a given ice loading history and an Earth model to calculate changes in sea level (i.e. a change in surface loading due to a change in the depth of the ocean) with time. The ocean load was computed using the ice loading histories W12, ICE-5G and ICE-6G C in combination with a three-layer Earth model and the output is a time-and space-variable load that can be applied to our laterally varying flat Earth model. We use an Earth structure that is representative of our 1-D average models with a lithospheric thickness of 96 km, upper-mantle viscosity of 5 × 10 20 Pa s, and lower-mantle viscosity of 1 × 10 22 Pa s. We acknowledge the inconsistencies inherent in this approach in that the ocean load is computed using a 1-D Earth model that may have different average upper-mantle viscosity and lithospheric thickness values to some of the models used in this study. However, we consider the impact of this to be small as there is, at most, ±0.7 mm yr −1 difference to present-day uplift rates when not including ocean loading at all; nevertheless, we choose to include ocean loading with the intention of making the model as realistic as possible. We keep the ocean loading the same for each ice model so that any differences in results may be attributed to differences in Earth properties. For the spatially uniform ice loading history, we do not include any ocean loading since it is an idealized loading history and would not produce a realistic sea-level change when modelled with a global GIA model.

R E S U LT S
In order to determine the effect of introducing a varying lithospheric thickness in experiment 1 (elastic case), we examine the differences between the 1-D and 3-D model output in terms of the spatial gradient of predicted present-day uplift rates. The spatial gradient is simply the derivative of the present-day uplift rate field and we take the scalar magnitude of the gradient (i.e. it is always positive) since the direction of the slope is not of interest. We calculate the spatial gradient over the 100 km resolution area of interest only. Differences in present-day uplift rates are relatively small (± 3 mm yr −1 , Fig. 3c) and the sign of the difference does not yield useful information. For example, in the Siple Coast the 3-D model predicts greater uplift at the coast but also more subsidence in the interior of West Antarctica, in other words, the 1-D model under-predicts the magnitude of the response compared with the 3-D model (Figs 3a and b). Differencing the deformation rates (1-D minus 3-D) shows both a positive and negative difference (Fig. 3c), masking the fact that the 3-D model produces a higher peak-to-peak difference in uplift rate (i.e. higher peaks and lower troughs). Calculating the spatial gradient of uplift rate for the 1-D and 3-D models and differencing them indicates how the amplitude and wavelength of deformation varies between the two models (Fig. 3d). A higher amplitude and shorter wavelength response would be expected from a thinner lithosphere compared with a thicker lithosphere (Wolstencroft et al. 2015). This can be observed in Fig. 4 from the profile of uplift rate and gradient of uplift rate across the Antarctic Peninsula, where the lithospheric thickness in the 3-D model is thinner than that of the 1-D average. The uplift rate predicted by the 3-D model (orange solid line in Fig. 4) has a higher amplitude and shorter wavelength (by one grid cell) than the 1-D model (green solid line in Fig. 4). This means that the gradient of uplift in the 3-D model will be steeper around the peak of the rebound, as indicated by the blue colour in the inset, but it tails off more quickly than in the 1-D model resulting in the 1-D gradient being steeper at the periphery (indicated by red in the inset). This gives a characteristic pattern of a white bulls eye (where the gradients are the same at the tip of the peak), surrounded by blue where the 3-D gradient is higher (negative gradient difference), with red at the periphery (Fig. 3d).
In East Antarctica, where the 1-D averaged lithospheric thickness is thinner than the 3-D model, and the present-day uplift rate gradients are steeper in the 1-D model output, the gradient difference is positive and shown as red, with the same characteristic white at the peak of the uplift/subsidence centres (Fig. 3d). Results for experiment 1 in Sections 3.1-3.3 are shown in the same format as Fig. 3(d)-as differences in uplift rate gradient between 1-D and 3-D models for our high resolution region of interest.  Fig. 5 shows the difference in uplift rate between these two models. The upper-and lower-mantle viscosities are kept the same for all models (5 × 10 20 and 1 × 10 22 Pa s, respectively). This plot can help us to understand what effect the ice loading history has on the results.

Effect of ice loading history
Each ice loading history results in different localized spatial patterns of present-day uplift rate gradient reflecting the spatial variability of ice loading or unloading between the models. Differences in the present-day uplift rate gradients of the 3-D and 1-D models, whether negative or positive, are focussed around the margins of centres (or 'peaks') of uplift or subsidence. This is because the lithospheric thickness in the 3-D model is thinner/thicker than the 1-D average and hence produces higher/lower amplitude and steeper/shallower gradients than the 1-D model. When comparing peaks of uplift rate gradient between the 1-D and 3-D model (for the same ice history) they have different amplitudes but the gradients at the crest of the peaks will often be the same or very similar, as explained previously, resulting in a small area of white at the centre of the region of uplift/subsidence (Fig. 5). For example, the ICE-6G C ice model (Figs 5c and g) shows a prominent blue bull'seye located near the Siple Coast related to a large unloading event. The unloading results in steeper uplift gradients and a higher peak amplitude in the 3-D case compared with the 1-D case, but in the centre, that is, at the peak itself, the gradient is the same (white on the figures). This effect can also be observed with the spatially uniform loading history (Figs 5d and h), where the periphery of the ice sheet shows the most sensitivity to variations in lithospheric thickness (i.e. largest differences in predicted present-day uplift rate gradients), and the interior shows little difference between the 1-D averaged model and the 3-D model.
Despite the localized differences in spatial pattern, all combinations of ice loading history and LAB model tested here yield the same first-order result across most of West Antarctica-use of a 1-D averaged lithospheric thickness results in lower magnitude gradients (lower amplitude and longer wavelength) of present-day uplift rate compared with the 3-D case, and hence predominantly negative differences across West Antarctica in Fig. 5. Any positive (red) differences in West Antarctica result from the longer wavelength deformation predicted by the 1-D model resulting in steeper gradients than the 3-D model at the periphery of the rebound. This result is insensitive to the ice model used, although the actual spatial patterns shown in Fig. 5 do depend on the ice loading history since the biggest differences in gradients when comparing a uniform lithospheric thickness to a laterally varying lithospheric thickness mostly occur around the margins of loading/unloading centres. The ice loading histories used in this study neglect any changes in icesheet thickness over the past few thousand years, such as those observed in the Antarctic Peninsula (Nield et al. 2012) and Siple Coast (Catania et al. 2012). Including these late Holocene changes would have the effect of changing the pattern of localized differences providing the underlying mantle viscosity was sufficiently low to respond on a timescale of ∼2000 yr or less.

Effect of LAB model
The choice of LAB model used to define spatial variations in lithospheric thickness has the potential to influence the results. The An et al. (2015a) model has a higher resolution than the Priestley & McKenzie (2013) model and therefore contains more spatial variability in the LAB depths. The bottom row of Figs 5 and 6 show the difference in uplift rate between the two LAB models, for the different ice loading models and upper-mantle viscosities, respectively. The impact of the LAB model in isolation can most clearly be observed in Fig. 5l-the model that uses the uniform loading history-because there are no spatial variations in ice loading that can amplify signals. The differences in Fig. 5l directly reflect the differences between the two LAB models (Fig. 1), with the greatest effects being in the Northern Antarctic Peninsula, where An et al. (2015a) identify a region of anomalously thick lithosphere, and in Coats Land (Fig. 2) where the boundary between East and West Antarctica is defined differently for each model. The differences peak at ±3.5 mm yr −1 for this latter region when using the W12 ice loading history (Fig. 5i) because large loading changes across this region during the past 5 ka (Whitehouse et al. 2012a) amplify the signal. All other ice loading/mantle viscosity combinations result in differences of ±1.5 mm yr −1 or less.
We can draw several conclusions from these observations. First, the results are more dependent on the ice loading history used than the choice of LAB model. Second, we do not gain significant extra information by using a higher resolution LAB model that resolves    the West having thinner-than-1-D-average lithosphere, as indicated by the dashed-dotted lines in Figs 5 and 6. This demarcation coincides with the regions where the amplitude of gradients of uplift rates for the 3-D model are higher (in West Antarctica) or lower (in East Antarctica) than the 1-D model and it is clearly the feature that has the most impact on gradients of uplift rates.

Effect of upper-mantle viscosity
Upper-mantle viscosity exerts a strong control on mantle relaxation times and hence uplift rates. To test whether our results are dependent on the underlying upper-mantle viscosity, we calculated the difference in present-day uplift rate gradients using four uppermantle viscosities, for both the LAB models in experiment 1, using just the W12 ice loading history (Fig. 6).
Comparing the results, we see similar patterns of gradient differences for the weaker upper-mantle viscosities (5 × 10 19 and 1 × 10 20 Pa s, Figs 6a and b, e and f) and the stronger upper-mantle viscosities (5 × 10 20 and 1 × 10 21 , Figs 6c and d, g and h), although the two sets of patterns are different from each other. The two sets of patterns reflect sensitivity to different periods in the deglacial history of the W12 ice model (Whitehouse et al. 2012a). The models with stronger mantle viscosities and slower relaxation time (Figs 6c  and d, g and h) are still rebounding in response to ice thinning in the western Ross Sea between 10 and 5 ka, whereas rebound in the lower viscosity models (Figs 6a and b, e and f) is dominated by the response to late Holocene ice thinning along the Siple Coast and Southern Antarctic Peninsula, as indicated by the blue areas in Figs 6(a and b) and (e and f). Fig. 6 demonstrates that the spatial variability in the gradient differences is dependent on both the ice loading history and the upper-mantle viscosity. Localized differences aside, for all viscosities we observe the same result of higher amplitude and shorter wavelength deformation in West Antarctica for the 3-D model (blue in the figures) supporting the hypothesis that the lithospheric thickness controls the wavelength of the signal captured in the modelling.

Effect of power-law rheology in the lithosphere
Modelling the lithosphere using a power-law rheology means that there is the potential for it to deform viscously, depending on the input temperature used to derive creep parameters and the stress from the ice loading. We compare results using power-law rheology (experiment 2) and input temperatures from An et al. (2015a; Fig. 7) with results from the equivalent experiment 1b model that has a spatially variable elastic lithosphere (Fig. 8); the two models have the same laterally varying lithospheric thickness but different rheology (see also Table 2). The upper-mantle viscosity (5 × 10 20 Pa s) and   ice loading history (W12) are the same for both models. Modelling the lithosphere with a power-law rheology has the effect of reducing the local effective elastic thickness (cf. Kuchar & Milne 2015) so we expect the power-law lithosphere (experiment 2) to behave as if it were thinner than the elastic lithosphere (experiment 1). In Fig. 8(a), we plot the difference in uplift rate gradient as elastic minus powerlaw so that the colour scale can be compared with the earlier plots of 1-D minus 3-D. The effective viscosities for elements that lie in the lithosphere are also calculated, following the methods described in van der Wal et al. (2015), and shown in Fig. 7 along with the temperatures from the An et al. (2015a) model that were used to derive the creep parameters.
The patterns of gradient difference show in Fig. 8(a) are unlike the previous results. Around the Weddell Sea (Fig. 2), there is a dark blue region indicating higher amplitude deformation in the powerlaw model compared with the elastic model in experiment 1, which may be related to the relatively low viscosity in the lithosphere at 70-90 km depth (around 1 × 10 22 Pa s, see the red circle in Fig. 7d) compared with the elastic lithosphere case (1 × 10 49 Pa s). Since the viscosity is dependent on the stress induced from ice load changes, the low viscosity in this region may also be associated with late ice loading changes defined within the W12 model. In fact, viscosity in this region is up to an order of magnitude lower (1 × 10 21 Pa s) during the load changes between 15 and 5 ka. Along the Siple Coast the large (blue) difference observed in the previous plots of Downloaded from https://academic.oup.com/gji/article-abstract/214/2/811/4980921 by University of Durham user on 11 June 2018 1-D versus 3-D is no longer present. This may be related to the fact that in this region, the seismic data indicate that there is a slab of relatively cold material at a depth of 50-70 km, resulting in a relatively high viscosity and therefore a very similar response to the case with the elastic lithosphere in experiment 1.
The profiles of present-day uplift rate and uplift rate gradient shown in Fig. 8(b) demonstrate that in the experiment that uses power-law rheology, the peaks have a higher amplitude and shorter wavelength than in the elastic lithosphere experiment. For the 50-70 km depth layer, the viscosity within the lower lithosphere under West Antarctica is around 1 × 10 20-21 Pa s, meaning it will deform viscously on glacial timescales of tens of thousands of years. This means that when using a power-law rheology to model the lithosphere only the upper 50 km will behave elastically over the timescales of interest (cf. Kuchar & Milne 2015).

Implications for future GIA models
In this study, we have shown that irrespective of deglacial history and sublithospheric mantle viscosity, the use of a spatially variable elastic lithospheric thickness in a GIA model of Antarctica results in higher gradients of predicted present-day uplift rates (i.e. higher amplitude and shorter wavelength) in West Antarctica compared with a uniform elastic lithospheric thickness that is simply the average of the former. We have made this comparison, first of all, to isolate the effect of introducing variable lithospheric thickness from any other factors that perturb predictions of uplift rates, and second, because many global GIA models use a 1-D Earth model derived from globally averaged parameters. The mean lithospheric thickness over the GIA model domain of both models of seismically derived LAB depth used in experiment 1 is 90 km, similar to values used in studies of global GIA [80-90 km, Mitrovica & Forte (2004) and Peltier et al. (2015), respectively]. Our results indicate that global 1-D GIA models with a ∼90 km lithospheric thickness would predict lower amplitude and longer wavelength uplift rates across West Antarctica than would be predicted with a more realistic, spatially variable lithosphere. This means that uplift rates, and hence geoid changes, would be smoothed out over a wider area potentially leading to an inaccurate GIA correction for GRACE data. A 1-D model with a lithospheric thickness representative of the average of West Antarctica (70 km) produces a closer match to results from the 3-D model than the 1-D Antarctic average lithosphere (90 km), apart from in regions where the lithosphere is even thinner (e.g. Southern Antarctic Peninsula, 50 km, Fig. 1c). This suggests that modelling East and West Antarctica with a separate 1-D Earth model is an important first step in improving GIA models of Antarctica.
Furthermore, modelling the lithosphere with a power-law rheology has the effect of reducing the thickness of the GIA lithosphere (i.e. the portion of the lithosphere acting elastically on GIA timescales) compared with the elastic case because the viscosity prescribed by the power-law rheology in the deeper parts of the lithosphere will be low enough to permit viscous behaviour over glacial timescales. By comparing results from experiment 1 and experiment 2, we have shown that using these two different definitions of the lithosphere leads to differences in gradients of presentday uplift rates despite input parameters (i.e. seismically derived LAB depth and seismically derived temperatures) ultimately coming from the same source. Using a power-law rheology provides a more consistent way of modelling GIA over multiple timescales because material properties determine the viscosity depending on timescale and this would, for example, allow relaxation of the lower lithosphere over multiple glacial cycles.
It is, therefore, important to consider both how the lithosphere is defined and how thickness variations are accounted for in the next generation of 3-D GIA models. As a minimum, East and West Antarctica should be considered separately in terms of Earth structure as both seismically derived LAB models considered here show a clear East-West divide in lithospheric thickness. We have shown that a model with higher resolution spatial variability in lithospheric thickness makes little difference to our results, however, representing lithospheric thickness variations within West Antarctica will become more important as ice loading histories evolve to contain greater spatial detail and include changes in ice thickness over the past few thousand years. Including a laterally varying lithospheric thickness would provide an improvement over current 1-D GIA models and should be considered to ensure more accurate predictions of uplift rate and ultimately a more accurate GIA correction for GRACE data. This is particularly pertinent for the dynamic region of West Antarctica that is currently experiencing a large amount of ice-mass loss (Rignot et al. 2014).

Implications for interpretation of observations of GIA
Geodetic observations of bedrock deformation provide useful data with which to constrain models of GIA. Consideration of laterally varying Earth structure may result in a better fit between model predictions and observations in some areas. For example, Wolstencroft et al. (2015) could not fit the spatial pattern of GPS-observed uplift in the southern Antarctic Peninsula with a 1-D Earth structure having tested several variations on recent ice loading history. It is possible that the strong spatial gradient in uplift revealed by differencing GPS rates recorded at sites on the east and west of the Antarctic Peninsula could be explained with the introduction of a thinner lithosphere in this region (e.g. 50-70 km, Fig. 1), which would be able to capture shorter wavelength differences in uplift, as we have shown. However, before such a comparison is made, Late Holocene ice loading changes (e.g. Nield et al. 2012;Nield et al. 2016) must be incorporated into current deglacial models.
Future observations of GIA should aim to be positioned in locations that would help to constrain 3-D Earth structure. In particular, increasing the density of GPS networks across West Antarctica would provide additional constraints for determining lithospheric thickness because shorter wavelength solid Earth deformation could be observed. For example, Nield et al. (2014) were able to more tightly constrain lithospheric thickness in the northern Antarctic Peninsula by using observations from the dense LARISSA network. Furthermore, measurements along the boundary between East and West Antarctica would provide useful information in delimiting this transition in Earth structure for the purposes of GIA models. Additional measurements of horizontal deformation could be instrumental in constraining lateral variations in Earth structure in this region.

Implications for ice-sheet models
We have demonstrated that the areas most affected by the inclusion of a spatially variable lithospheric thickness lie around the margins of ice loading changes, including (for most combinations of ice history and Earth model tested) the Amundsen Sea sector and the Siple Coast (locations shown in Fig. 2). This has important implications Downloaded from https://academic.oup.com/gji/article-abstract/214/2/811/4980921 by University of Durham user on 11 June 2018 for ice dynamics in marine-grounded areas that lie on a reverse slope bed, for example, West Antarctica. Grounding line dynamics control ice-sheet stability and evidence shows that a reverse slope bed can reduce ice-sheet stability because, as the grounding line retreats into deeper water, ice flux across the grounding line will increase, potentially leading to net ice loss and hence further retreat (e.g. Schoof 2007).
Studies of Antarctic ice loss that make use of a coupled icesheet-sea-level model have shown that bedrock uplift has a stabilizing effect on marine-grounded ice due to reducing the slope of a reverse bed, resulting in less ice loss from Antarctica (Gomez et al. 2010;Gomez et al. 2013). Including a spatially variable lithospheric thickness would increase the stabilizing effect of bedrock uplift on the marine-grounded sector of the ice sheet in West Antarctica compared with a 1-D averaged model because, as we have shown, the thinner lithosphere results in higher amplitude uplift in the interior, thereby reducing the slope of the reverse bed further. This has been demonstrated by Gomez et al. (2015) and Pollard et al. (2017) who show that a 1-D Earth model with a 50 km lithospheric thickness and low mantle viscosity results in increased stabilization over a 1-D model with thicker lithosphere and higher mantle viscosity. Furthermore, Gomez et al. (2018) showed that a coupled ice-sheet-sea-level model with a 3-D Earth structure (laterally varying lithospheric thickness and upper-mantle viscosity) results in significant regional differences in ice-sheet thickness when compared with results using a 1-D Earth structure. In particular, their model predicts thicker ice and less retreat of the grounding line over the last deglaciation at the periphery of the Ross Sea region (Fig. 2) where the lithosphere is thinner, and upper-mantle viscosity is lower, than their 1-D average model. Including 3-D Earth structure in GIA models and ice dynamic models is, therefore, necessary for determining the dynamics of past ice-sheet change and accurately assessing the current and future state of the West Antarctic Ice Sheet.

Limitations
Model resolution is an important consideration for any GIA model. Here, we restricted the spatial resolution to 100 × 100 km elements in the area of interest, purely for computational efficiency. We tested the effect of running a higher resolution model, increasing the mesh resolution to 50 km in the area of interest. While the output is smoother, the 50 km resolution model did not reveal any additional features that are not captured by the 100 km mesh and considering the extra computation time, we conclude that the coarser resolution is satisfactory for the experiments performed in this study.
In the finite element model, material properties are considered compressible in the computation of deformation, but the effect this has on buoyancy forces is not considered. The model also neglects self-gravitation, that is, changes in gravitational potential caused by deformation, which is a feature of most spherical models. However, Schotman et al. (2008) state that when using a flat-earth model, the lack of sphericity partly cancels the lack of self-gravitation. Furthermore, since we are looking at differences between models, any errors arising due to the neglect of such features will effectively be cancelled out.

C O N C L U S I O N S
We have presented the results of two experiments that seek to investigate the impact of including lateral variations in lithospheric thickness when modelling the solid Earth response to surface loading across West Antarctica. The first experiment used estimates for the depth of the LAB derived from seismic studies to model the lithosphere as an elastic layer, an approach taken in many GIA studies. We have compared results from 3-D models (varying lithospheric thickness) and equivalent 1-D models (uniform lithospheric thickness that is the average of the 3-D model). For all combinations of ice history, LAB model and underlying upper-mantle viscosity tested, we find that the use of a 1-D averaged lithospheric thickness results in lower gradients (i.e. lower amplitude and longer wavelength) of uplift rate compared with the use of a spatially variable (thinner in West Antarctica) lithospheric thickness. This means that the present-day uplift rate is smoothed over a wider area in the 1-D model and the magnitude of peaks and troughs of deformation is smaller. This has important implications for ice-sheet modelling studies as steeper spatial gradients of uplift may promote stability in marine-grounded regions of West Antarctica.
The biggest difference in results between the two different seismically derived LAB models used is in the Northern Antarctic Peninsula and at the boundary between East and West Antarctica, partly due to the An et al. (2015a) model having higher resolution and a greater level of detail. The most important feature of these LAB models is the delineation of where the lithosphere is thinner than average in West Antarctica, which is a stable feature across different seismic models, although the location of this boundary is important as it can affect uplift rates in this area. Within West Antarctica the localized patterns of differences in uplift rate gradient are sensitive to the choice of ice loading history, with largest differences focussed around centres of loading or unloading. The choice of underlying mantle viscosity also plays a role because the viscosity defines the relaxation time of the mantle, which in turn determines which regions will still be deforming in response to past ice-sheet change. Including a laterally variable lithospheric thickness within West Antarctica will become even more important once ice loading histories incorporate changes from the past few thousand years.
The second experiment in this study investigated the difference between two methods of defining the lithosphere in GIA modelling. We compared the elastic lithosphere in experiment 1 with the use of power-law rheology in experiment 2, which defines viscosity based on material parameters and loading changes, and hence implicitly defines the lithosphere based on whether the viscosity is high enough to behave elastically over the timescale in question. Our results demonstrate that using a power-law rheology produces higher amplitude peaks of deformation than using a 3-D elasticonly lithosphere because in the power-law case, the thickness of the portion of the lithosphere that behaves elastically is reduced. Defining the lithosphere in this way could provide a more robust model of GIA since the thickness of the lithosphere is less rigidly defined than in the elastic (i.e. very high viscosity) case and relaxation in the lower lithosphere could be important when modelling several glacial cycles (Kaufmann et al. 2005).
Future GIA models should seek to include a spatially varying lithospheric thickness, or at the very least to represent thinner/thicker lithosphere in West/East Antarctica; we find that inclusion of this transition has a first order effect on the predicted pattern of present-day deformation. Regional 1-D GIA models should ensure that the local lithospheric thickness is adequately represented rather than using an average of a wider Antarctic domain. Furthermore, including a spatially variable lithosphere could lead to a better fit to GPS-observed uplift rates, especially in regions where a thinner lithosphere might be necessary to capture shore wavelength Downloaded from https://academic.oup.com/gji/article-abstract/214/2/811/4980921 by University of Durham user on 11 June 2018 signals. This, in turn, could improve GIA models in West Antarctica where the uncertainty is large, although lateral variations in mantle viscosity and better constraints on ice history would also be required to provide an improved correction for GRACE data.

A C K N O W L E D G E M E N T S
We thank Rebekka Steffen and an anonymous reviewer for their constructive comments that helped to improve the manuscript. GAN and JPOD are supported by NERC grant NE/L006294/1. PLW is a recipient of an NERC Independent Research Fellowship (NE/K009958/1). This research is a contribution to the SCAR SERCE program. All figures have been produced using the GMT software package (Wessel & Smith 1998).