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

Grace A Nield; Pippa L Whitehouse; Wouter van der Wal; Bas Blank; John Paul O'Donnell; Graham W Stuart

doi:10.1093/gji/ggy158

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

Nield, Grace A; Whitehouse, Pippa L; van der Wal, Wouter; Blank, Bas; O'Donnell, John Paul; Stuart, Graham W 2018-08-01 00:00:00 Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Geophys. J. Int. (2018) 214, 811–824 doi: 10.1093/gji/ggy158 Advance Access publication 2018 April 19 GJI Geodynamics and tectonics The impact of lateral variations in lithospheric thickness on glacial isostatic adjustment in West Antarctica 1 1 2 2 Grace A. Nield, Pippa L. Whitehouse, Wouter van der Wal, Bas Blank, John 3 3 Paul O’Donnell and Graham W. Stuart Department of Geography, Durham University, Durham DH1 3LE, UK. E-mail: grace.a.nield@durham.ac.uk Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS, Delft, The Netherlands School of Earth and Environment, University of Leeds, Leeds LS2 9JT, UK Accepted 2018 April 19. Received 2018 April 17; in original form 2018 January 30 SUMMARY 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. Key words: Creep and deformation; Satellite geodesy; Antarctica; Dynamics of lithosphere and mantle; Rheology: crust and lithosphere; Rheology: mantle. due to large uncertainties that persist in both the ice-sheet history 1. INTRODUCTION since the Last Glacial Maximum (LGM) and the Earth structure in The process of Glacial Isostatic Adjustment (GIA) in Antarctica is this region. This has a direct impact on estimates of ice-mass loss well-studied (e.g. Whitehouse et al. 2012b;A et al. 2013; Argus derived from satellite gravimetry (e.g. the Gravity Recovery and et al. 2014) but GIA models continue to predict remarkably dif- Climate Experiment, GRACE) since Antarctic GIA is a signiﬁcant ferent present-day deformation rates (Mart´ ın-Espanol ˜ et al. 2016) The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 812 G. A. Nield et al. component of the total gravitational signal and must be removed to 2005; Nield et al. 2014; Wolstencroft et al. 2014). On the basis of yield estimates for ice-mass balance (King et al. 2012). wave speed variations, seismic studies can distinguish between ther- Traditionally, global and Antarctic-wide models of GIA have mal conduction and convection regimes in the upper mantle. The used a 1-D approximation of the Earth structure consisting of an conductive domain deﬁnes the tectonic plate. However, the elastic elastic lithosphere underlain by a linear viscoelastic upper and lower thickness varies as a function of timescale of surface loading and is mantle, where properties vary only radially (e.g. Peltier 1974; Milne typically thinner than the seismic lithospheric thickness. & Mitrovica 1996; Kendall et al. 2005). In reality, the structure of The studies and methods described above have used linear vis- the Earth is far more complex and models that reﬂect lateral as well coelastic rheology to model GIA. However, the use of power-law as vertical variations in Earth properties are needed to provide more rheology is becoming increasingly common (Wu 1999; Barnhoorn accurate predictions of present-day GIA-related deformation and et al. 2011; van der Wal et al. 2013; van der Wal et al. 2015). van geoid changes, both in Antarctica (A et al. 2013;van derWal et al. der Wal et al. (2015) used seismic velocity anomalies (Grand 2002) 2015;Sasgen et al. 2017) and elsewhere, for example, Greenland and geothermal heat ﬂux estimates (Shapiro & Ritzwoller 2004) (Khan et al. 2016). Including 3-D structure in GIA models is par- for Antarctica to infer mantle temperatures that were used to derive ticularly pertinent for Antarctica as this continent is considered to creep parameters for input to a power-law rheology. By deﬁning spa- consist of two distinct regions in terms of Earth structure: a thick tially varying creep parameters, the GIA model included laterally cratonic lithosphere and a high-viscosity uppermost mantle in the varying Earth structure. For this approach, the lithospheric thick- East, and thinner lithosphere and lower viscosity uppermost mantle ness is implicitly deﬁned by the creep parameters, rheological model in the West (Morelli & Danesi 2004). Modelling East and West and some threshold viscosity above which it can be considered to Antarctica with a 1-D Earth model as described above, therefore, behave elastically as described above. has the potential to produce incorrect estimates of the present- Modelling advances in the past few decades (Wu & Johnston day GIA signal in one or both of these regions. For example, A 1998; Latychev et al. 2005b;A et al. 2013;van derWal et al. 2013) et al. (2013) compared deformation rates predicted by a 3-D model have eased the computational burden of 3-D GIA modelling and incorporating laterally varying lithospheric thickness and mantle detailed data sets are now available that can be used to deﬁne lateral viscosity with a model that is the 1-D average of the 3-D proﬁle Earth structure (Ritsema et al. 2011;Heeszel et al. 2016). Hence, and found mismatches at Global Positioning System (GPS) loca- there are an increased number of studies incorporating 3-D Earth tions in Antarctica. Furthermore, capturing variability in the Earth structure into GIA models with both linear and nonlinear rheologies. structure within West Antarctica is important because regional 1-D Several approaches can be used to infer 3-D mantle viscosity (Ivins GIA studies have indicated differences in the Earth structure across & Sammis 1995;Kaufmann et al. 2005) and lithospheric thickness the region (Nield et al. 2014; Wolstencroft et al. 2015; Nield et al. for input to GIA models, with the latter being the focus of this 2016;Zhao et al. 2017). study. A seismically derived lithosphere–asthenosphere boundary This study focusses on how lateral variations in lithospheric thick- (LAB) depth is sometimes used to infer laterally varying GIA litho- ness impact predictions made by GIA models. The lithospheric spheric thickness with linear viscoelastic rheology, after scaling to thickness can be deﬁned by various criteria, such as a change in the account for differences between a seismically derived deﬁnition of method of heat transfer (Martinec & Wolf 2005), seismic anisotropy the lithosphere and the mechanical deﬁnition used in GIA studies. or resistivity (Eaton et al. 2009). In GIA modelling, the lithosphere For example, Kaufmann et al. (2000) reduced a seismically derived is deﬁned on the basis of mechanical properties and is considered LAB depth by a factor of two for their GIA modelling, and Khan to be a part of the crust and upper mantle that behaves elastically et al. (2016) used an adjustment factor to scale the LAB depth pub- on timescales of glacial cycles (tens of thousands of years; Mar- lished by Priestley & McKenzie (2013). However, it is not clear tinec & Wolf 2005; Watts et al. 2013; Kuchar & Milne 2015). The whether a lithosphere deﬁned by seismic properties can be con- lithosphere can be modelled with either a purely elastic rheology, verted to a lithosphere deﬁned by mechanical properties through a that is, it has no viscous component (e.g. Argus et al. 2014), or as a scaling factor. Seismic properties could have a different dependence viscoelastic material with sufﬁciently high viscosity that it does not on temperature and composition than mechanical properties. One way to circumvent this issue is to use temperatures derived from relax in response to surface loading on timescales of a glacial cycle (e.g. Kendall et al. 2005; Whitehouse et al. 2012b; Kuchar & Milne seismic velocity perturbations as input to a power-law rheology, 2015), thereby behaving elastically. Studies have also combined which eliminates the need to explicitly deﬁne lithospheric thick- these approaches, for example, Kaufmann et al. (2005) modelled a ness (van der Wal et al. 2013). In this approach, assumptions are 100 km thick lithosphere composed of a 30 km purely elastic layer made in converting seismic velocity anomalies into temperature and overlying a 70 km viscoelastic layer with a viscosity of 1 × 10 Pa viscosity, and the lithosphere is deﬁned implicitly by the effective s, which is approximately the limit of what could be considered viscosity. So although temperatures ultimately come from the same elastic over GIA timescales [e.g. 1 × 10 Pa s (Sasgen et al. 2017) seismic source as the LAB depths, no new assumptions are required –1 × 10 Pa s (Khan et al. 2016)]. Kuchar & Milne (2015)in- other than those made for converting seismic velocities to viscosity. vestigated the effect of depth-dependent viscosity in the lithosphere Previous studies investigating 3-D Earth structure in GIA mod- on relative sea-level predictions using a radially varying (i.e. 1-D) els of Antarctica have focussed on the effect of lateral variations Earth model and found that predictions made using a lithosphere in mantle viscosity (e.g. Kaufmann et al. 2005) or a combination with viscosity structure were similar to predictions made using a of laterally varying lithospheric thickness and upper-mantle viscos- purely elastic, but much thinner, lithosphere. ity (A et al. 2013; van der Wal et al. 2015) on present-day uplift To some extent, the apparent thickness of the lithosphere de- rates. Studies that isolate the effect of including lateral variations pends on the timescales of surface loading. Over long timescales in lithospheric thickness in models of GIA exist for regions in the (∼1 Myr), viscous relaxation in the lower lithosphere means that the northern hemisphere (Kaufmann et al. 2000; Zhong et al. 2003; lithosphere seems to behave as a relatively thin elastic layer (Watts Latychev et al. 2005a; Whitehouse et al. 2006;Steffen et al. 2014) et al. 2013). However, over GIA timescales (∼100 kyr), the litho- but currently no such study exists for Antarctica. sphere seems to behave as a thicker elastic layer (Martinec & Wolf Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 813 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 deﬁning a laterally varying lithospheric thickness mentioned above. The ﬁrst method (experiment 1) uses a scaled seismically inferred LAB depth to de- termine 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 differ- ences 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 litho- spheric 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 differ- ences 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 tem- perature as input to a power-law relationship, thereby assigning a viscosity structure to the lithosphere. In this method, the lithospheric thickness is implicitly deﬁned as the depth at which the resulting viscosity is small enough so that signiﬁcant deformation takes place Figure 1. Adjusted elastic lithospheric thickness derived from (a) Priestley during a glacial cycle. In order to determine the effect of includ- & McKenzie (2013)and (c)An et al. (2015a) LAB depths. Each colour ing viscosity in the lithosphere through a power-law rheology on represents a separate 20 km thick model layer with the lithospheric thickness being the upper bound of the colour, for example, orange denotes a LAB gradients of present-day uplift rate, we compare results using the depth/lithospheric thickness of 90 km. (b) and (d) show where the 3-D power-law rheology to those from the ﬁrst method which assumes lithosphere is thinner or thicker than the 90 km 1-D average. The regular that the lithosphere is elastic. For both methods, we use a ﬁnite mesh of 100 × 100 km is bounded by locations (−3000 km, −2000 km) element model with 20 km thick layers, meaning that the variable and (500 km, 2500 km), with an irregular mesh outside of this region. lithospheric thickness is captured in 20 km ‘steps’ between element locations. Due to large uncertainties in both Earth structure and at the edge of the model domain, for computational efﬁciency, and ice history, we do not attempt to ﬁt any observational data such as the domain has an overall width of 60 000 km. The extremely large GPS-observed uplift. model domain is required to ensure that boundary effects are negli- gible in the area of interest (Steffen et al. 2006). We do not model any ice loading changes outside Antarctica as the impact on uplift 2 METHODS AND DATA rates in Antarctica would be negligible (Whitehouse et al. 2012b) but we do include ocean loading. The model consists of 22 depth 2.1 Model and geometry layers representing the Earth’s surface to the core–mantle boundary (Table 1). A 30 km purely elastic crustal layer [the same thickness We use a 3-D ﬂat-earth ﬁnite element model constructed with the used by Kaufmann et al. (2005)] is underlain by eleven 20 km thick ABAQUS software package (Hibbitt et al. 2016) to compute the layers to 250 km depth to capture a higher resolution in the litho- solid Earth deformation in response to a changing surface load. sphere and uppermost mantle. Below this, layers are thicker (i.e. The validity of using this approach to model the Earth’s response lower resolution with depth) since surface deformation will be less to changes in ice-sheet loading has been shown previously by Wu sensitive to the details of mantle rheology below 250 km depth (Lau & Johnston (1998). This method has been used in many studies to et al. 2016). The buoyancy force is accounted for by applying Win- model GIA in regions such as Fennoscandia (Kaufmann et al. 2000; kler foundations to layer boundaries where a density contrast occurs Steffen et al. 2006), Antarctica (Kaufmann et al. 2005) and Iceland (Wu 2004). To ensure a direct comparability between the models, (Auriac et al. 2013), and has the advantage of computational efﬁ- the same mesh is used for all the experiments. ciency over a spherical global model. The ﬂat-earth ﬁnite element approach has been shown to be accurate when computing defor- mation within the ice margin for ice loads with comparable size to 2.2 Earth models and data the Laurentide Ice Sheet (Wu & Johnston 1998), which makes it applicable to the Antarctic Ice Sheet with its smaller lateral extent. The compressible elastic material properties for each layer described The mesh consists of eight-node brick elements with reduced above are listed in Table 1. The elastic and density structures of integration. The surface geometry of the mesh consists of a the Earth are derived from the preliminary reference Earth model 3500 × 4500 km area of interest embedded in a larger model (Dziewonski & Anderson 1981). For each element below the upper- domain. The area of interest represents West Antarctica and has most elastic layer, creep parameters are assigned on an element-by- elements of 100 × 100 km (elements are shown in Fig. 1). Outside element basis. The geometry of the mesh means that the lithospheric this region, the element size increases towards the periphery of the thickness (experiment 1) or the viscosity (experiment 2) can vary in model, from 550 km in East Antarctica to approximately 5000 km 20 km steps between the adjacent element locations. Latychev et al. Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 814 G. A. Nield et al. Table 1. Model layers and Earth parameters. Top of layer Top of layer Layer thickness Young’s modulus − 3 Layer radius (km) depth (km) (km) Density (kg m ) (GPa) Poisson’s ratio Rheology Lithosphere 6371 0 30 3196 173.9 0.28 Elastic Lithosphere or 6341 30 20 3379 173.9 0.28 Elastic/Power-law UM 6321 50 20 3377 173.9 0.28 lithosphere, or 6301 70 20 3375 173.3 0.28 linear viscoelastic 6281 90 20 3373 172.7 0.28 if UM 6261 110 20 3370 171.6 0.28 6241 130 20 3368 170.6 0.28 6221 150 20 3366 170.0 0.28 6201 170 20 3364 169.3 0.28 6181 190 20 3362 179.5 0.29 6161 210 20 3436 194.6 0.29 6141 230 20 3448 200.8 0.30 6121 250 80 3478 212.6 0.30 UM 6041 330 70 3525 224.4 0.30 Linear 5971 400 136 3812 277.2 0.29 viscoelastic- 5835 536 134 3978 377.8 0.28 variable LM 5701 670 251 4482 459.4 0.27 Linear viscoelastic 5450 921 250 4630 484.2 0.28 1 × 10 Pa s 5200 1171 430 4825 509.0 0.28 4770 1601 430 5036 570.1 0.29 4340 2031 430 5264 636.9 0.29 3910 2461 430 5464 704.5 0.30 (2005a) demonstrated that the differences in uplift rate over previ- Kaufmann et al. 2005; Kendall et al. 2005; Whitehouse et al. 2012b; ously glaciated regions in the northern hemisphere peak at 2 mm Ivins et al. 2013; Wolstencroft et al. 2015) and for the timescales −1 yr for a jump of 150 km between continental and oceanic litho- we are interested in the lithosphere is generally considered to be 24 25 spheric thickness, so we conclude that the effect of a 20 km jump elastic at viscosities above 1 × 10 –1 × 10 Pa s (Kaufmann et al. on the uplift rate is likely to be small. Our approach to deﬁning vari- 2005;Steffen et al. 2006; Barnhoorn et al. 2011;Khan et al. 2016). able lithospheric thickness allows us to use the same mesh for all Throughout the rest of this paper, we simply refer to this type of experiments, thus ensuring that the results are directly comparable. modelled lithosphere as the ‘elastic lithosphere’. For each model The sublithospheric upper mantle (down to a depth of 670 km) is a in experiment 1, we compare the laterally varying model with an linear viscoelastic layer with uniform viscosity and several different equivalent 1-D model in which the lithospheric thickness is simply upper-mantle viscosities are tested to determine dependence of the the average of the 3-D lithospheric thickness. results on the underlying viscosity (see Table2). Properties for the The second experiment uses a power-law rheology; this comple- lower mantle are the same for all models (Table 1). mentary approach allows us to investigate the differences in the The thickness of the lithosphere is deﬁned differently in ex- two methods used to deﬁne the lithospheric thickness in GIA mod- periments 1 and 2, as detailed in Table 2. The ﬁrst experiment elling. Mantle temperatures (An et al. 2015a) are used to determine considers an elastic lithosphere with spatial variability deﬁned by diffusion and dislocation creep parameters following the methods two different models of seismically derived LAB depth (Priestley described by Hirth & Kohlstedt (2003)and vander Wal et al. (2013; & McKenzie 2013;An et al. 2015a), as described in the follow- 2015). The reader is referred to these papers for a detailed descrip- ing sections. Seismically derived LAB depths tend to be thicker tion of the method. We limit the power-law rheology to the same than those inferred from GIA studies. For example, Steffen et al. horizontal and vertical domain as deﬁned by the An et al. (2015a) (2014) compare GIA-inferred lithospheric thicknesses in the Baltic LAB depths for two reasons: (1) so that the results of experiment 2 Sea with three LAB depth models and ﬁnd a consistently thinner can be directly compared with the results of experiment 1, and (2) lithosphere by 30–80 km. We, therefore, uniformly reduce the LAB so that the upper-mantle viscosity remains laterally uniform and, depths from the Priestley & McKenzie (2013)and An et al. (2015a) therefore, the effect of a spatially variable lithospheric thickness is models so that the thicknesses are more representative of a GIA isolated from all other parameters. 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 & 2.2.1 Priestley & McKenzie (2013) McKenzie (2013) because it is an Antarctic-speciﬁc model derived Priestley & McKenzie (2013) published a model of global seismic using many additional seismic stations. For those elements repre- velocities from the surface wave tomography which they used to senting the lithosphere, the viscosity is set to 1 × 10 Pa s to mimic derive mantle temperatures and lithospheric thickness on a 2 grid elastic behaviour on glacial timescales, apart from the uppermost with a resolution of 250 km horizontally and 50 km vertically. The 30 km layer which is modelled as purely elastic. This approach lithospheric thickness is deﬁned by the change in heat transfer from of modelling an elastic lithosphere using a viscoelastic rheology conduction to convection. We use this information in experiment 1a with high viscosity (including combinations of purely elastic and to deﬁne an elastic lithospheric thickness. The authors state that the viscoelastic rheology) is taken in many GIA studies (Peltier 2004; uncertainty on the lithospheric thickness is 20–30 km, and therefore Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 815 Table 2. Summary of the experiments and inputs used. Lithosphere Experiment deﬁnition Data used Ice models Upper-mantle viscosity (Pa s) Results 1a Elastic Priestley & McKenzie (2013) W12 (all viscosities); ICE-5G, 5 × 10 (all ice models); Comparison between 1-D and 19 20 21 LAB depths (adjusted) ICE-6G C, Uniform loading 5 × 10 ,1 × 10 ,1 × 10 3-D (for 5 × 10 Pa s) (W12 only) 1b Elastic An et al. (2015a) LAB depths W12 (all viscosities); ICE-5G, 5 × 10 (all ice models); Comparison between 1-D and 19 20 21 (adjusted) ICE-6G C, Uniform loading 5 × 10 ,1 × 10 ,1 × 10 3-D (for 5 × 10 Pa s) (W12 only) 2 Power-law Domain deﬁned by An et al. W12 5 × 10 Comparison between elastic (2015a) LAB depths, An et al. 3-D (1b) and power-law 3-D (2015b) temperatures used as (same ice model and input to power-law rheology upper-mantle viscosity) we reduce the LAB depths by 20 km to reﬂect the fact that a GIA- The LAB depths mapped onto the model elements is shown in inferred elastic lithosphere is typically thinner than a lithosphere Fig. 1(c), again with Fig. 1(d) showing where the 3-D model has based on the change in heat transfer method (Martinec & Wolf thicker or thinner lithosphere than the 1-D average [90 km, calcu- 2005; Priestley & McKenzie 2013). GIA-inferred LAB thicknesses lated over the Antarctic Plate which is the spatial limit of the inferred are less than the LAB depths in the Priestley & McKenzie (2013) LABdepthsinthe An et al. (2015a) model]. Similar to Priestley & model for regions such as Iceland [less than 50 per cent, 15–40 km McKenzie (2013), the lithosphere under East Antarctica is thicker (Auriac et al. 2013) compared with ∼95 km] and Fennoscandia than the 1-D average for the An et al. (2015a) model. The location [around 50–70 per cent, 93–110 km (Zhao et al. 2012) compared of the boundary between East and West Antarctica is similar in both with 120–200 km], and we use the uncertainty bound to reduce the seismic models, indicating that the uncertainty on the location of LAB depths to 70–90 per cent of the modelled values so that our the boundary is small. Some isolated regions of anomalously thick results are conservative. Fig. 1(a) shows the adjusted lithospheric lithosphere are also present in the Northern Antarctic Peninsula, thicknesses mapped onto the ABAQUS layers (Table 1), that is, at which the authors attribute to a remnant subducted slab from the each location on the mesh, the adjusted LAB depth is rounded to former subduction zone in this region. 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 2.3 Ice loading a thicker than average lithosphere. 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 2.2.2 An et al. (2015a,b) et al. 2013; Argus et al. 2014). Given this uncertainty, one of the The second model used in this study is from An et al. (2015a), who aims of this study is to investigate whether differences between 1- infer temperatures below Antarctica from the 3-D seismic velocity D and 3-D models are independent of the assumed ice history. To model AN1-S (An et al. 2015b), which has a horizontal resolution test this, we use several different deglacial models, or ice loading that increases from ∼120 km in the crust to ∼500 km at a depth of histories, that have quite different spatial patterns and magnitude 120 km and a vertical resolution of ∼25 to 150 km depth followed of loading changes, which, when applied to a speciﬁc Earth model, by ∼50 to 250 km depth. We use this model in two ways. First, the give different patterns of deformation. We compare gradients of seismically derived LAB, which is deﬁned by the depth where the present-day uplift rates between 1-D and 3-D Earth models using adiabat crosses the 1330 C geotherm, is used to deﬁne lithospheric the same ice history, revealing differences that may be directly at- thickness for the elastic lithosphere in experiment 1b. The uncer- tributed to the introduction of a varying lithospheric thickness, and tainty on the temperature is reported to be ±150 C, equivalent to then qualitatively compare the results from different ice models. ±15–30 km for the LAB depth, so we reduce the LAB depth by Three ice loading scenarios are used in the modelling: W12 15 km to be representative of GIA-elastic thicknesses, as explained (Whitehouse et al. 2012a), ICE-5G (Peltier 2004) and its successor in Section 1, and to be consistent with the scaling of the Priestley ICE-6G C (Argus et al. 2014; Peltier et al. 2015). The W12 ice & McKenzie (2013) model. Second, we use the temperatures di- loading model was derived using a glaciologically consistent nu- rectly as input to the power-law rheology in experiment 2. In this merical ice-sheet model that was tuned to provide the best possible experiment, we consider the 3-D spatial domain that deﬁnes the ﬁt to constraints of ice thickness change, whereas ICE-5G and ICE- lithosphere in experiment 1, but within this domain we use a power 6G C have been tuned to ﬁt ice thickness change constraints and law instead of elastic rheology. Results show the comparison of the GPS-observed uplift rates without satisfying ice-sheet physics. Fur- two 3-D models, thereby highlighting the differences due to rheo- thermore, in an attempt to fully isolate differences associated with logical deﬁnitions (see Table 2 for a summary of the models). In the introduction of a laterally varying lithospheric thickness from their model An et al. (2015a) do not infer temperatures for depths those caused by spatial variations in ice loading, we also construct shallower than 55 km. Therefore, when using the temperatures in an idealistic, spatially uniform loading history. In this scenario, the our model we specify a second elastic layer between 30 and 50 km amount of ice thickness change since the LGM is spatially uniform depth to bridge the gap between our uppermost elastic layer and the over the grounded area of the present-day ice sheet (as shown in temperature inputs. Table3 and Fig. 2) with a somewhat arbitrary 700 m of total ice loss Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 816 G. A. Nield et al. Table 3. Ice thickness change for the spatially uniform ice loading history for the West Antarctic Ice Sheet (WAIS) and the East Antarctic Ice Sheet (EAIS). Ice thickness change (m) Time period (ka) WAIS EAIS 20–15 −200 50 15–10 −300 60 10–5 −150 30 5–0 −50 10 Total: LGM to Present −700 150 Figure 3. Present-day uplift rates for the a) 1-D and b) 3-D models based on An et al. (2015a) LAB depths (experiment 1b), using the W12 ice loading history and upper-mantle viscosity 5 × 10 Pa s; c) difference in present- day uplift rates (1-D minus 3-D); d) difference in spatial gradient of uplift rate between 1-D and 3-D model (1-D minus 3-D) for the high resolution region of interest only—blue areas show where the 3-D model predicts higher amplitude and shorter wavelength deformation compared with the Figure 2. Regions of ice thickness change for the uniform loading his- 1-D model. tory for West Antarctica (blue) and East Antarctica (green). Key locations mentioned in the text also labelled. lithospheric thickness values to some of the models used in this study. However, we consider the impact of this to be small as there for West Antarctica and 150 m ice-sheet growth for East Antarc- −1 is, at most, ±0.7 mm yr difference to present-day uplift rates tica, applied over four time periods. These ice thickness changes when not including ocean loading at all; nevertheless, we choose approximate the mean ice loading changes in the W12 ice loading to include ocean loading with the intention of making the model as history [compare with ﬁg. 7 of Whitehouse et al. (2012a)]. realistic as possible. We keep the ocean loading the same for each ice model so that any differences in results may be attributed to dif- ferences in Earth properties. For the spatially uniform ice loading 2.4 Ocean loading history, we do not include any ocean loading since it is an idealized The model approach we have used in this study does not solve loading history and would not produce a realistic sea-level change the sea-level equation (Farrell & Clark 1976) and cannot compute when modelled with a global GIA model. variable sea level with time. We, therefore, take the approach of applying an ocean load that has been derived using a global, spher- ically symmetric GIA model (Mitrovica & Milne 2003; Kendall 3 RESULTS et al. 2005; Mitrovica et al. 2005). The GIA model uses a given ice loading history and an Earth model to calculate changes in In order to determine the effect of introducing a varying lithospheric sea level (i.e. a change in surface loading due to a change in the thickness in experiment 1 (elastic case), we examine the differences depth of the ocean) with time. The ocean load was computed using between the 1-D and 3-D model output in terms of the spatial the ice loading histories W12, ICE-5G and ICE-6G C in combina- gradient of predicted present-day uplift rates. The spatial gradient tion with a three-layer Earth model and the output is a time- and is simply the derivative of the present-day uplift rate ﬁeld and we space-variable load that can be applied to our laterally varying ﬂat take the scalar magnitude of the gradient (i.e. it is always positive) Earth model. We use an Earth structure that is representative of since the direction of the slope is not of interest. We calculate the our 1-D average models with a lithospheric thickness of 96 km, spatial gradient over the 100 km resolution area of interest only. upper-mantle viscosity of 5 × 10 Pa s, and lower-mantle viscosity Differences in present-day uplift rates are relatively small 22 −1 of 1 × 10 Pa s. We acknowledge the inconsistencies inherent in (± 3mmyr , Fig. 3c) and the sign of the difference does not this approach in that the ocean load is computed using a 1-D Earth yield useful information. For example, in the Siple Coast the 3-D model that may have different average upper-mantle viscosity and model predicts greater uplift at the coast but also more subsidence Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 817 in the interior of West Antarctica, in other words, the 1-D model ICE-6G C ice model (Figs 5c and g) shows a prominent blue bull’s- under-predicts the magnitude of the response compared with the eye located near the Siple Coast related to a large unloading event. 3-D model (Figs 3a and b). Differencing the deformation rates (1-D The unloading results in steeper uplift gradients and a higher peak minus 3-D) shows both a positive and negative difference (Fig. 3c), amplitude in the 3-D case compared with the 1-D case, but in the masking the fact that the 3-D model produces a higher peak-to-peak centre, that is, at the peak itself, the gradient is the same (white difference in uplift rate (i.e. higher peaks and lower troughs). Calcu- on the ﬁgures). This effect can also be observed with the spatially lating the spatial gradient of uplift rate for the 1-D and 3-D models uniform loading history (Figs 5d and h), where the periphery of and differencing them indicates how the amplitude and wavelength the ice sheet shows the most sensitivity to variations in lithospheric of deformation varies between the two models (Fig. 3d). A higher thickness (i.e. largest differences in predicted present-day uplift rate amplitude and shorter wavelength response would be expected from gradients), and the interior shows little difference between the 1-D a thinner lithosphere compared with a thicker lithosphere (Wolsten- averaged model and the 3-D model. croft et al. 2015). This canbeobservedinFig. 4 from the proﬁle of Despite the localized differences in spatial pattern, all combi- uplift rate and gradient of uplift rate across the Antarctic Peninsula, nations of ice loading history and LAB model tested here yield where the lithospheric thickness in the 3-D model is thinner than the same ﬁrst-order result across most of West Antarctica—use of that of the 1-D average. The uplift rate predicted by the 3-D model a 1-D averaged lithospheric thickness results in lower magnitude (orange solid line in Fig. 4) has a higher amplitude and shorter gradients (lower amplitude and longer wavelength) of present-day wavelength (by one grid cell) than the 1-D model (green solid line uplift rate compared with the 3-D case, and hence predominantly in Fig. 4). This means that the gradient of uplift in the 3-D model negative differences across West Antarctica in Fig. 5. Any positive will be steeper around the peak of the rebound, as indicated by the (red) differences in West Antarctica result from the longer wave- blue colour in the inset, but it tails off more quickly than in the 1-D length deformation predicted by the 1-D model resulting in steeper model resulting in the 1-D gradient being steeper at the periphery gradients than the 3-D model at the periphery of the rebound. This (indicated by red in the inset). This gives a characteristic pattern result is insensitive to the ice model used, although the actual spa- of a white bulls eye (where the gradients are the same at the tip tial patterns shown in Fig. 5 do depend on the ice loading history of the peak), surrounded by blue where the 3-D gradient is higher since the biggest differences in gradients when comparing a uniform (negative gradient difference), with red at the periphery (Fig. 3d). lithospheric thickness to a laterally varying lithospheric thickness In East Antarctica, where the 1-D averaged lithospheric thickness mostly occur around the margins of loading/unloading centres. The is thinner than the 3-D model, and the present-day uplift rate gra- ice loading histories used in this study neglect any changes in ice- dients are steeper in the 1-D model output, the gradient difference sheet thickness over the past few thousand years, such as those is positive and shown as red, with the same characteristic white observed in the Antarctic Peninsula (Nield et al. 2012) and Siple at the peak of the uplift/subsidence centres (Fig. 3d). Results for Coast (Catania et al. 2012). Including these late Holocene changes experiment 1 in Sections 3.1–3.3 are shown in the same format as would have the effect of changing the pattern of localized differ- Fig. 3(d)—as differences in uplift rate gradient between 1-D and ences providing the underlying mantle viscosity was sufﬁciently 3-D models for our high resolution region of interest. low to respond on a timescale of ∼2000 yr or less. 3.2 Effect of LAB model 3.1 Effect of ice loading history The choice of LAB model used to deﬁne spatial variations in litho- Fig. 5 shows the difference in spatial gradient of the present-day spheric thickness has the potential to inﬂuence the results. The An uplift rate when comparing the 3-D and 1-D models for the four et al. (2015a) model has a higher resolution than the Priestley & different ice loading histories used in this study. Results are shown McKenzie (2013) model and therefore contains more spatial vari- for both models of lithospheric thickness used in experiment 1 – ability in the LAB depths. The bottom row of Figs 5 and 6 show Priestley & McKenzie (2013)and An et al. (2015a); the bottom the difference in uplift rate between the two LAB models, for the row in Fig. 5 shows the difference in uplift rate between these two different ice loading models and upper-mantle viscosities, respec- models. The upper- and lower-mantle viscosities are kept the same tively. The impact of the LAB model in isolation can most clearly 20 22 for all models (5 × 10 and 1 × 10 Pa s, respectively). This plot be observed in Fig. 5l—the model that uses the uniform loading can help us to understand what effect the ice loading history has on history—because there are no spatial variations in ice loading that the results. can amplify signals. The differences in Fig. 5l directly reﬂect the Each ice loading history results in different localized spatial pat- differences between the two LAB models (Fig. 1), with the greatest terns of present-day uplift rate gradient reﬂecting the spatial vari- effects being in the Northern Antarctic Peninsula, where An et al. ability of ice loading or unloading between the models. Differences (2015a) identify a region of anomalously thick lithosphere, and in in the present-day uplift rate gradients of the 3-D and 1-D mod- Coats Land (Fig. 2) where the boundary between East and West els, whether negative or positive, are focussed around the margins Antarctica is deﬁned differently for each model. The differences −1 of centres (or ‘peaks’) of uplift or subsidence. This is because peak at ±3.5 mm yr for this latter region when using the W12 ice the lithospheric thickness in the 3-D model is thinner/thicker than loading history (Fig. 5i) because large loading changes across this the 1-D average and hence produces higher/lower amplitude and region during the past 5 ka (Whitehouse et al. 2012a) amplify the steeper/shallower gradients than the 1-D model. When comparing signal. All other ice loading/mantle viscosity combinations result −1 peaks of uplift rate gradient between the 1-D and 3-D model (for in differences of ±1.5 mm yr or less. the same ice history) they have different amplitudes but the gradi- We can draw several conclusions from these observations. First, ents at the crest of the peaks will often be the same or very similar, the results are more dependent on the ice loading history used than as explained previously, resulting in a small area of white at the the choice of LAB model. Second, we do not gain signiﬁcant extra centre of the region of uplift/subsidence (Fig. 5). For example, the information by using a higher resolution LAB model that resolves Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 818 G. A. Nield et al. Figure 4. Uplift rate (left-hand axis) for the 1-D (solid green) and 3-D (solid orange) models along the proﬁle shown in the inset. Also shown is the gradient of uplift rate (right hand axis) along the proﬁle for the 1-D (dashed green) and 3-D (dashed orange) models, with shading according to the difference in gradient shown in the inset (1-D minus 3-D; same as Fig. 3d). Black dashed line indicates the difference in gradient shown in the inset plot. Figure 5. Difference in spatial gradient of uplift rate between 1-D and 3-D models (1-D minus 3-D) for ice loading histories (from left to right) W12 (a, e), ICE-5G (b, f), ICE-6G C (c, g) and the uniform loading history (d, h), and for the two different LAB models, Priestley & McKenzie (2013;top row) andAn et al. (2015a; middle row). All models have an upper-mantle viscosity of 5 × 10 Pa s. The dashed-dotted black line delineates where the 3-D lithosphere is thinner or thicker than in the 1-D case, as shown in Figs 1(b) and (d). Panels (i)-(l) show the difference in uplift rate between the 3-D LAB models (Priestley & McKenzie (2013) minus An et al. (2015a)). smaller scale variations in lithospheric thickness, even if we increase both seismically inferred LAB models show a clear East–West di- the GIA model horizontal resolution to 50 km (Section 4.4). Finally, vide, with the East having thicker-than-1-D-average lithosphere and Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 819 Figure 6. Difference in spatial gradient of uplift rate between 1-D and 3-D models (1-D minus 3-D) for different values of upper-mantle viscosity (from left to right), for the two different LAB models, Priestley & McKenzie (2013;top row) andAn et al. (2015a; middle row), and using the W12 ice history. The dashed-dotted black line delineates where the 3-D lithosphere is thinner or thicker than in the 1-D case, as shown in Figs 1(b) and (d). Panels (i)–(l) show the difference in uplift rate between the 3-D models for the two different LAB models [Priestley & McKenzie (2013) minus An et al. (2015a)]. the West having thinner-than-1-D-average lithosphere, as indicated the western Ross Sea between 10 and 5 ka, whereas rebound in the by the dashed-dotted lines in Figs 5 and 6. This demarcation co- lower viscosity models (Figs 6a and b, e and f) is dominated by incides with the regions where the amplitude of gradients of uplift the response to late Holocene ice thinning along the Siple Coast rates for the 3-D model are higher (in West Antarctica) or lower (in and Southern Antarctic Peninsula, as indicated by the blue areas in East Antarctica) than the 1-D model and it is clearly the feature that Figs 6(a and b) and (e and f). has the most impact on gradients of uplift rates. 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 viscosi- ties we observe the same result of higher amplitude and shorter 3.3 Effect of upper-mantle viscosity wavelength deformation in West Antarctica for the 3-D model (blue Upper-mantle viscosity exerts a strong control on mantle relaxation in the ﬁgures) supporting the hypothesis that the lithospheric thick- times and hence uplift rates. To test whether our results are de- ness controls the wavelength of the signal captured in the modelling. pendent on the underlying upper-mantle viscosity, we calculated the difference in present-day uplift rate gradients using four upper- mantle viscosities, for both the LAB models in experiment 1, using 3.4 Effect of power-law rheology in the lithosphere just the W12 ice loading history (Fig. 6). Comparing the results, we see similar patterns of gradient dif- Modelling the lithosphere using a power-law rheology means that ferences for the weaker upper-mantle viscosities (5 × 10 and there is the potential for it to deform viscously, depending on the 1 × 10 Pa s, Figs 6a and b, e and f) and the stronger upper-mantle input temperature used to derive creep parameters and the stress 20 21 viscosities (5 × 10 and 1 × 10 , Figs 6c and d, g and h), although from the ice loading. We compare results using power-law rheology the two sets of patterns are different from each other. The two sets of (experiment 2) and input temperatures from An et al. (2015a; Fig. 7) patterns reﬂect sensitivity to different periods in the deglacial his- with results from the equivalent experiment 1b model that has a spa- tory of the W12 ice model (Whitehouse et al. 2012a). The models tially variable elastic lithosphere (Fig. 8); the two models have the with stronger mantle viscosities and slower relaxation time (Figs 6c same laterally varying lithospheric thickness but different rheology and d, g and h) are still rebounding in response to ice thinning in (see also Table 2). The upper-mantle viscosity (5 × 10 Pa s) and Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 820 G. A. Nield et al. Figure 7. Top row: temperatures from the An et al. (2015b) model averaged over the ﬁnite element model layers. Bottom row: effective viscosity at the present-day for the same model layers as the top row, calculated following the methods detailed in van der Wal et al. (2015). Red circle in panel (d) shows low viscosity lithosphere mentioned in the text. Elements below the spatially variable lithospheric thickness from An et al. (2015a) are greyed out (cf. Fig. 1c). Figure 8. (a) Difference in spatial gradient between the 3-D elastic-only case (experiment 1b) and the 3-D power-law case (experiment 2; elastic-only case minus power-law case), for the W12 ice loading history with upper-mantle viscosity of 5 × 10 Pa s. (b) Proﬁle of uplift rate for the elastic (green solid) and power-law (orange solid) cases and the gradient of each (dashed lines, right hand axis) along the proﬁle shown in (a). ice loading history (W12) are the same for both models. Modelling The patterns of gradient difference show in Fig. 8(a) are unlike the lithosphere with a power-law rheology has the effect of reducing the previous results. Around the Weddell Sea (Fig. 2), there is a dark the local effective elastic thickness (cf. Kuchar & Milne 2015)sowe blue region indicating higher amplitude deformation in the power- expect the power-law lithosphere (experiment 2) to behave as if it law model compared with the elastic model in experiment 1, which were thinner than the elastic lithosphere (experiment 1). In Fig. 8(a), may be related to the relatively low viscosity in the lithosphere at we plot the difference in uplift rate gradient as elastic minus power- 70–90 km depth (around 1 × 10 Pa s, see the red circle in Fig. 7d) law so that the colour scale can be compared with the earlier plots compared with the elastic lithosphere case (1 × 10 Pa s). Since the of 1-D minus 3-D. The effective viscosities for elements that lie in viscosity is dependent on the stress induced from ice load changes, the lithosphere are also calculated, following the methods described the low viscosity in this region may also be associated with late ice in van der Wal et al. (2015), and shown in Fig. 7 along with the loading changes deﬁned within the W12 model. In fact, viscosity temperatures from the An et al. (2015a) model that were used to in this region is up to an order of magnitude lower (1 × 10 Pa derive the creep parameters. 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/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 821 1-D versus 3-D is no longer present. This may be related to the because material properties determine the viscosity depending on fact that in this region, the seismic data indicate that there is a slab timescale and this would, for example, allow relaxation of the lower of relatively cold material at a depth of 50–70 km, resulting in a lithosphere over multiple glacial cycles. relatively high viscosity and therefore a very similar response to the It is, therefore, important to consider both how the lithosphere case with the elastic lithosphere in experiment 1. is deﬁned and how thickness variations are accounted for in the The proﬁles of present-day uplift rate and uplift rate gradient next generation of 3-D GIA models. As a minimum, East and West shown in Fig. 8(b) demonstrate that in the experiment that uses Antarctica should be considered separately in terms of Earth struc- power-law rheology, the peaks have a higher amplitude and shorter ture as both seismically derived LAB models considered here show wavelength than in the elastic lithosphere experiment. For the 50– a clear East-West divide in lithospheric thickness. We have shown 70 km depth layer, the viscosity within the lower lithosphere under that a model with higher resolution spatial variability in lithospheric 20–21 West Antarctica is around 1 × 10 Pa s, meaning it will de- thickness makes little difference to our results, however, represent- form viscously on glacial timescales of tens of thousands of years. ing lithospheric thickness variations within West Antarctica will This means that when using a power-law rheology to model the become more important as ice loading histories evolve to contain lithosphere only the upper 50 km will behave elastically over the greater spatial detail and include changes in ice thickness over the timescales of interest (cf. Kuchar & Milne 2015). 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 predic- tions of uplift rate and ultimately a more accurate GIA correction for 4 DISCUSSION GRACE data. This is particularly pertinent for the dynamic region of West Antarctica that is currently experiencing a large amount of 4.1 Implications for future GIA models ice-mass loss (Rignot et al. 2014). 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 4.2 Implications for interpretation of observations of GIA in higher gradients of predicted present-day uplift rates (i.e. higher Geodetic observations of bedrock deformation provide useful data amplitude and shorter wavelength) in West Antarctica compared with which to constrain models of GIA. Consideration of laterally with a uniform elastic lithospheric thickness that is simply the av- varying Earth structure may result in a better ﬁt between model erage of the former. We have made this comparison, ﬁrst of all, to predictions and observations in some areas. For example, Wolsten- isolate the effect of introducing variable lithospheric thickness from croft et al. (2015) could not ﬁt the spatial pattern of GPS-observed any other factors that perturb predictions of uplift rates, and second, uplift in the southern Antarctic Peninsula with a 1-D Earth struc- because many global GIA models use a 1-D Earth model derived ture having tested several variations on recent ice loading history. from globally averaged parameters. The mean lithospheric thickness It is possible that the strong spatial gradient in uplift revealed by over the GIA model domain of both models of seismically derived differencing GPS rates recorded at sites on the east and west of the LAB depth used in experiment 1 is 90 km, similar to values used Antarctic Peninsula could be explained with the introduction of a in studies of global GIA [80–90 km, Mitrovica & Forte (2004)and thinner lithosphere in this region (e.g. 50–70 km, Fig. 1), which Peltier et al. (2015), respectively]. Our results indicate that global would be able to capture shorter wavelength differences in uplift, as 1-D GIA models with a ∼90 km lithospheric thickness would pre- we have shown. However, before such a comparison is made, Late dict lower amplitude and longer wavelength uplift rates across West Holocene ice loading changes (e.g. Nield et al. 2012; Nield et al. Antarctica than would be predicted with a more realistic, spatially 2016) must be incorporated into current deglacial models. variable lithosphere. This means that uplift rates, and hence geoid Future observations of GIA should aim to be positioned in loca- changes, would be smoothed out over a wider area potentially lead- tions that would help to constrain 3-D Earth structure. In particular, ing to an inaccurate GIA correction for GRACE data. A 1-D model increasing the density of GPS networks across West Antarctica with a lithospheric thickness representative of the average of West would provide additional constraints for determining lithospheric Antarctica (70 km) produces a closer match to results from the 3-D thickness because shorter wavelength solid Earth deformation could model than the 1-D Antarctic average lithosphere (90 km), apart be observed. For example, Nield et al. (2014)wereabletomore from in regions where the lithosphere is even thinner (e.g. Southern tightly constrain lithospheric thickness in the northern Antarctic Antarctic Peninsula, 50 km, Fig. 1c). This suggests that modelling Peninsula by using observations from the dense LARISSA network. East and West Antarctica with a separate 1-D Earth model is an Furthermore, measurements along the boundary between East and important ﬁrst step in improving GIA models of Antarctica. West Antarctica would provide useful information in delimiting Furthermore, modelling the lithosphere with a power-law rhe- this transition in Earth structure for the purposes of GIA models. ology has the effect of reducing the thickness of the GIA litho- Additional measurements of horizontal deformation could be in- sphere (i.e. the portion of the lithosphere acting elastically on GIA strumental in constraining lateral variations in Earth structure in timescales) compared with the elastic case because the viscosity this region. 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 4.3 Implications for ice-sheet models experiment 2, we have shown that using these two different deﬁni- tions of the lithosphere leads to differences in gradients of present- We have demonstrated that the areas most affected by the inclusion day uplift rates despite input parameters (i.e. seismically derived of a spatially variable lithospheric thickness lie around the margins LAB depth and seismically derived temperatures) ultimately com- of ice loading changes, including (for most combinations of ice his- ing from the same source. Using a power-law rheology provides tory and Earth model tested) the Amundsen Sea sector and the Siple a more consistent way of modelling GIA over multiple timescales Coast (locations shown in Fig. 2). This has important implications Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 822 G. A. Nield et al. for ice dynamics in marine-grounded areas that lie on a reverse thickness when modelling the solid Earth response to surface load- slope bed, for example, West Antarctica. Grounding line dynamics ing across West Antarctica. The ﬁrst experiment used estimates for control ice-sheet stability and evidence shows that a reverse slope the depth of the LAB derived from seismic studies to model the bed can reduce ice-sheet stability because, as the grounding line lithosphere as an elastic layer, an approach taken in many GIA retreats into deeper water, ice ﬂux across the grounding line will studies. We have compared results from 3-D models (varying litho- increase, potentially leading to net ice loss and hence further retreat spheric thickness) and equivalent 1-D models (uniform lithospheric (e.g. Schoof 2007). thickness that is the average of the 3-D model). For all combinations Studies of Antarctic ice loss that make use of a coupled ice- of ice history, LAB model and underlying upper-mantle viscosity sheet–sea-level model have shown that bedrock uplift has a stabi- tested, we ﬁnd that the use of a 1-D averaged lithospheric thickness lizing effect on marine-grounded ice due to reducing the slope of a results in lower gradients (i.e. lower amplitude and longer wave- reverse bed, resulting in less ice loss from Antarctica (Gomez et al. length) of uplift rate compared with the use of a spatially variable 2010; Gomez et al. 2013). Including a spatially variable lithospheric (thinner in West Antarctica) lithospheric thickness. This means that thickness would increase the stabilizing effect of bedrock uplift on the present-day uplift rate is smoothed over a wider area in the 1-D the marine-grounded sector of the ice sheet in West Antarctica com- model and the magnitude of peaks and troughs of deformation is pared with a 1-D averaged model because, as we have shown, the smaller. This has important implications for ice-sheet modelling thinner lithosphere results in higher amplitude uplift in the inte- studies as steeper spatial gradients of uplift may promote stability rior, thereby reducing the slope of the reverse bed further. This in marine-grounded regions of West Antarctica. has been demonstrated by Gomez et al. (2015) and Pollard et al. The biggest difference in results between the two different seis- (2017) who show that a 1-D Earth model with a 50 km lithospheric mically derived LAB models used is in the Northern Antarctic thickness and low mantle viscosity results in increased stabiliza- Peninsula and at the boundary between East and West Antarctica, tion over a 1-D model with thicker lithosphere and higher mantle partly due to the An et al. (2015a) model having higher resolution viscosity. Furthermore, Gomez et al. (2018) showed that a cou- and a greater level of detail. The most important feature of these pled ice-sheet–sea-level model with a 3-D Earth structure (laterally LAB models is the delineation of where the lithosphere is thinner varying lithospheric thickness and upper-mantle viscosity) results than average in West Antarctica, which is a stable feature across in signiﬁcant regional differences in ice-sheet thickness when com- different seismic models, although the location of this boundary pared with results using a 1-D Earth structure. In particular, their is important as it can affect uplift rates in this area. Within West model predicts thicker ice and less retreat of the grounding line Antarctica the localized patterns of differences in uplift rate gradi- over the last deglaciation at the periphery of the Ross Sea region ent are sensitive to the choice of ice loading history, with largest (Fig. 2) where the lithosphere is thinner, and upper-mantle viscosity differences focussed around centres of loading or unloading. The is lower, than their 1-D average model. Including 3-D Earth struc- choice of underlying mantle viscosity also plays a role because ture in GIA models and ice dynamic models is, therefore, necessary the viscosity deﬁnes the relaxation time of the mantle, which in for determining the dynamics of past ice-sheet change and accu- turn determines which regions will still be deforming in response rately assessing the current and future state of the West Antarctic to past ice-sheet change. Including a laterally variable lithospheric Ice Sheet. thickness within West Antarctica will become even more important once ice loading histories incorporate changes from the past few thousand years. 4.4 Limitations The second experiment in this study investigated the difference between two methods of deﬁning the lithosphere in GIA modelling. Model resolution is an important consideration for any GIA model. We compared the elastic lithosphere in experiment 1 with the use Here, we restricted the spatial resolution to 100 × 100 km elements of power-law rheology in experiment 2, which deﬁnes viscosity in the area of interest, purely for computational efﬁciency. We tested based on material parameters and loading changes, and hence im- the effect of running a higher resolution model, increasing the mesh plicitly deﬁnes the lithosphere based on whether the viscosity is resolution to 50 km in the area of interest. While the output is high enough to behave elastically over the timescale in question. smoother, the 50 km resolution model did not reveal any additional Our results demonstrate that using a power-law rheology produces features that are not captured by the 100 km mesh and considering higher amplitude peaks of deformation than using a 3-D elastic- the extra computation time, we conclude that the coarser resolution only lithosphere because in the power-law case, the thickness of the is satisfactory for the experiments performed in this study. portion of the lithosphere that behaves elastically is reduced. Deﬁn- In the ﬁnite element model, material properties are considered ing the lithosphere in this way could provide a more robust model compressible in the computation of deformation, but the effect this of GIA since the thickness of the lithosphere is less rigidly deﬁned has on buoyancy forces is not considered. The model also neglects than in the elastic (i.e. very high viscosity) case and relaxation in self-gravitation, that is, changes in gravitational potential caused by the lower lithosphere could be important when modelling several deformation, which is a feature of most spherical models. However, glacial cycles (Kaufmann et al. 2005). Schotman et al. (2008) state that when using a ﬂat-earth model, Future GIA models should seek to include a spatially vary- the lack of sphericity partly cancels the lack of self-gravitation. ing lithospheric thickness, or at the very least to represent thin- Furthermore, since we are looking at differences between models, ner/thicker lithosphere in West/East Antarctica; we ﬁnd that inclu- any errors arising due to the neglect of such features will effectively sion of this transition has a ﬁrst order effect on the predicted pattern be cancelled out. of present-day deformation. Regional 1-D GIA models should en- sure that the local lithospheric thickness is adequately represented rather than using an average of a wider Antarctic domain. Fur- 5 CONCLUSIONS thermore, including a spatially variable lithosphere could lead to a We have presented the results of two experiments that seek to in- better ﬁt to GPS-observed uplift rates, especially in regions where a vestigate the impact of including lateral variations in lithospheric thinner lithosphere might be necessary to capture shore wavelength Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 823 signals. This, in turn, could improve GIA models in West Antarctica Heeszel, D.S. et al., 2016. Upper mantle structure of central and West Antarc- tica from array analysis of Rayleigh wave phase velocities, J. geophys. where the uncertainty is large, although lateral variations in mantle Res., 121, 1758–1775. viscosity and better constraints on ice history would also be required Hibbitt, D., Karlsson, B. & Sorensen, P., 2016. Getting Started with ABAQUS, to provide an improved correction for GRACE data. Version (6.14), Hibbitt, Karlsson & Sorensen, Inc. Hirth, G. & Kohlstedt, D., 2003. Rheology of the upper mantle and the mantle wedge: a view from the experimentalists, in Inside the Subduction ACKNOWLEDGEMENTS Factory, ed. Eiler, J., pp. 83–105, American Geophysical Union. Ivins, E.R., James, T.S., Wahr, J., Schrama, E.J.O., Landerer, F.W. & Simon, We thank Rebekka Steffen and an anonymous reviewer for their K.M., 2013. Antarctic contribution to sea level rise observed by GRACE constructive comments that helped to improve the manuscript. with improved GIA correction, J. geophys. Res., 118, 3126–3141. GAN and JPOD are supported by NERC grant NE/L006294/1. Ivins, E.R. & Sammis, C.G., 1995. On lateral viscosity contrast in the mantle PLW is a recipient of an NERC Independent Research Fellow- and the rheology of low-frequency geodynamics, Geophys. J. Int., 123, ship (NE/K009958/1). This research is a contribution to the SCAR 305–322. SERCE program. All ﬁgures have been produced using the GMT Kaufmann, G., Wu, P. & Ivins, E.R., 2005. Lateral viscosity variations software package (Wessel & Smith 1998). beneath Antarctica and their implications on regional rebound motions and seismotectonics, J. Geodyn., 39, 165–181. Kaufmann, G., Wu, P. & Li, G., 2000. Glacial isostatic adjustment in Fennoscandia for a laterally heterogeneous earth, Geophys. J. Int., 143, REFERENCES 262–273. A, G, Wahr, J. & Zhong, S., 2013. Computations of the viscoelastic response Kendall, R.A., Mitrovica, J.X. & Milne, G.A., 2005. On post-glacial sea of a 3-D compressible Earth to surface loading: an application to Glacial level – II. Numerical formulation and comparative results on spherically Isostatic Adjustment in Antarctica and Canada, Geophys. J. Int., 192, symmetric models, Geophys. J. Int., 161, 679–706. 557–572. Khan, S.A. et al., 2016. Geodetic measurements reveal similarities between An, M. et al., 2015a. Temperature, lithosphere–asthenosphere boundary, and post–last glacial maximum and present-day mass loss from the Greenland heat ﬂux beneath the Antarctic Plate inferred from seismic velocities, J. ice sheet, Sci. Adv., 2, doi:10.1126/sciadv.1600931. geophys. Res., 120, 8720–8742. King, M.A., Bingham, R.J., Moore, P., Whitehouse, P.L., Bentley, M.J. & An, M. et al., 2015b. S-velocity model and inferred Moho topography Milne, G.A., 2012. Lower satellite-gravimetry estimates of Antarctic sea- beneath the antarctic plate from rayleigh waves, J. geophys. Res., 120, level contribution, Nature, 491, 586–589. 359–383. Kuchar, J. & Milne, G.A., 2015. The inﬂuence of viscosity structure in the Argus, D.F., Peltier, W.R., Drummond, R. & Moore, A.W., 2014. The Antarc- lithosphere on predictions from models of glacial isostatic adjustment, J. tica component of postglacial rebound model ICE-6G C (VM5a) based Geodyn., 86, 1–9. on GPS positioning, exposure age dating of ice thicknesses, and relative Latychev, K., Mitrovica, J.X., Tamisiea, M.E., Tromp, J. & Moucha, R., sea level histories, Geophys. J. Int., 198, 537–563 2005a. Inﬂuence of lithospheric thickness variations on 3-D crustal veloc- Auriac, A., Spaans, K.H., Sigmundsson, F., Hooper, A., Schmidt, P. & Lund, ities due to glacial isostatic adjustment, Geophys. Res. Lett., 32, L01304. B., 2013. Iceland rising: solid Earth response to ice retreat inferred from Latychev, K., Mitrovica, J.X., Tromp, J., Tamisiea, M.E., Komatitsch, D. & satellite radar interferometry and visocelastic modeling, J. geophys. Res., Christara, C.C., 2005b. Glacial isostatic adjustment on 3-D Earth models: 118, 1331–1344. a ﬁnite-volume formulation, Geophys. J. Int., 161, 421–444. Barnhoorn, A., van der Wal, W., Vermeersen, B.L.A. & Drury, M.R., 2011. Lau, H.C.P., Mitrovica, J.X., Austermann, J., Crawford, O., Al-Attar, D. & Lateral, radial, and temporal variations in upper mantle viscosity and Latychev, K., 2016. Inferences of mantle viscosity based on ice age data rheology under Scandinavia, Geochem. Geophys. Geosyst., 12, Q01007, sets: Radial structure, J. geophys. Res., 121, 6991–7012. doi:10.1029/2010GC003290. Mart´ ın-Espanol, ˜ A., King, M.A., Zammit-Mangion, A., Andrews, S.B., Briggs, R.D. & Tarasov, L., 2013. How to evaluate model-derived deglacia- Moore, P. & Bamber, J.L., 2016. An assessment of forward and inverse tion chronologies: a case study using Antarctica, Quat. Sci. Rev., 63, GIA solutions for Antarctica, J. geophys. Res., 121, 6947–6965. 109–127. Martinec, Z. & Wolf, D., 2005. Inverting the Fennoscandian relaxation- Catania, G., Hulbe, C., Conway, H., Scambos, T.A. & Raymond, C.F., 2012. time spectrum in terms of an axisymmetric viscosity distribution with a Variability in the mass ﬂux of the Ross ice streams, West Antarctica, over lithospheric root, J. Geodyn., 39, 143–163. the last millennium, J. Glaciol., 58, 741–752. Milne, G.A. & Mitrovica, J.X., 1996. Postglacial sea-level change on a Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference Earth rotating Earth: ﬁrst results from a gravitationally self-consistent sea-level model, Phys. Earth planet. Inter., 25, 297–356. equation, Geophys. J. Int., 126, F13–F20. Eaton, D.W., Darbyshire, F., Evans, R.L., Grutter, H., Jones, A.G. & Yuan, X., Mitrovica, J.X. & Forte, A.M., 2004. A new inference of mantle viscosity 2009. The elusive lithosphere–asthenosphere boundary (LAB) beneath based upon joint inversion of convection and glacial isostatic adjustment cratons, Lithos, 109, 1–22. data, Earth planet. Sci. Lett., 225, 177–189. Farrell, W.E. & Clark, J.A., 1976. On postglacial sea level, Geophys. J. R. Mitrovica, J.X. & Milne, G.A., 2003. On post-glacial sea level: I. General astr. Soc., 46, 647–667. theory, Geophys. J. Int., 154, 253–267. Gomez, N., Latychev, K. & Pollard, D., 2018. A coupled ice sheet-sea level Mitrovica, J.X., Wahr, J., Matsuyama, I. & Paulson, A., 2005. The rotational model incorporating 3D Earth structure: variations in Antarctica during stability of an ice-age earth, Geophys. J. Int., 161, 491–506. the last deglacial retreat, J. Clim., doi:10.1175/JCLI-D-17-0352.1. Morelli, A. & Danesi, S., 2004. Seismological imaging of the Antarctic Gomez, N., Mitrovica, J.X., Huybers, P. & Clark, P.U., 2010. Sea level as continental lithosphere: a review, Glob. Planet. Change, 42, 155–165. a stabilizing factor for marine-ice-sheet grounding lines, Nat. Geosci, 3, Nield, G.A. et al., 2014. Rapid bedrock uplift in the Antarctic Peninsula 850–853. explained by viscoelastic response to recent ice unloading, Earth planet. Gomez, N., Pollard, D. & Holland, D., 2015. Sea-level feedback lowers Sci. Lett., 397, 32–41. projections of future Antarctic Ice-Sheet mass loss, Nat. Commun., 6, Nield, G.A., Whitehouse, P.L., King, M.A. & Clarke, P.J., 2016. Glacial 8798, doi:10.1038/ncomms9798. isostatic adjustment in response to changing Late Holocene behaviour of Gomez, N., Pollard, D. & Mitrovica, J.X., 2013. A 3-D coupled ice sheet – ice streams on the Siple Coast, West Antarctica, Geophys. J. Int., 205, sea level model applied to Antarctica through the last 40 ky, Earth planet. 1–21. Sci. Lett., 384, 88–99. Nield, G.A., Whitehouse, P.L., King, M.A., Clarke, P.J. & Bentley, M.J., Grand, S.P., 2002. Mantle shear–wave tomography and the fate of subducted 2012. Increased ice loading in the Antarctic Peninsula since the 1850s slabs, Phil. Trans. R. Soc. Lond., A, 360, 2475–2491. Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 824 G. A. Nield et al. and its effect on glacial isostatic adjustment, Geophys. Res. Lett., 39, van der Wal, W., Whitehouse, P.L. & Schrama, E.J.O., 2015. Effect of GIA L17504. models with 3D composite mantle viscosity on GRACE mass balance Peltier, W.R., 1974. The impulse response of a Maxwell Earth, Rev. Geophys., estimates for Antarctica, Earth planet. Sci. Lett., 414, 134–143. 12, 649–669. Watts, A.B., Zhong, S.J. & Hunter, J., 2013. The behavior of the lithosphere Peltier, W.R., 2004. Global glacial isostasy and the surface of the ice-age on seismic to geologic timescales, Annu. Rev. Earth planet. Sci., 41, 443– Earth: the ICE-5G (VM2) Model and GRACE, Annu. Rev. Earth planet. 468. Sci., 32, 111–149. Wessel, P. & Smith, W.H.F., 1998. New, improved version of generic mapping Peltier, W.R., Argus, D.F. & Drummond, R., 2015. Space geodesy constrains tools released, EOS, Trans. Am. geophys. Un., 79, 579–579. ice age terminal deglaciation: the global ICE-6G C (VM5a) model, J. Whitehouse, P., Latychev, K., Milne, G.A., Mitrovica, J.X. & Kendall, R., geophys. Res., 120, 450–487. 2006. Impact of 3-D Earth structure on Fennoscandian glacial isostatic ad- Pollard, D., Gomez, N. & Deconto, R.M., 2017. Variations of the Antarctic justment: Implications for space-geodetic estimates of present-day crustal ice sheet in a coupled ice sheet-Earth-Sea level model: sensitivity to deformations, Geophys. Res. Lett., 33, doi:10.1029/2006GL026568. viscoelastic Earth properties, J. geophys. Res., 122, 2124–2138. Whitehouse, P.L., Bentley, M.J. & Le Brocq, A.M., 2012a. A deglacial Priestley, K. & McKenzie, D., 2013. The relationship between shear wave model for Antarctica: geological constraints and glaciological modelling velocity, temperature, attenuation and viscosity in the shallow part of the as a basis for a new model of Antarctic glacial isostatic adjustment, Quat. mantle, Earth planet. Sci. Lett., 381, 78–91. Sci. Rev., 32, 1–24. Rignot, E., Mouginot, J., Morlighem, M., Seroussi, H. & Scheuchl, B., 2014. Whitehouse, P.L., Bentley, M.J., Milne, G.A., King, M.A. & Thomas, I.D., Widespread, rapid grounding line retreat of Pine Island, Thwaites, Smith, 2012b. A new glacial isostatic adjustment model for Antarctica: calibrated and Kohler glaciers, West Antarctica, from 1992 to 2011, Geophys. Res. and tested using observations of relative sea-level change and present-day Lett., 41, 3502–3509. uplift rates, Geophys. J. Int., 190, 1464–1482. Ritsema, J., Deuss, A., van Heijst, H.J. & Woodhouse, J.H., 2011. S40RTS: Wolstencroft, M. et al., 2015. Uplift rates from a new high-density GPS a degree-40 shear-velocity model for the mantle from new Rayleigh wave network in Palmer Land indicate signiﬁcant late Holocene ice loss in the dispersion, teleseismic traveltime and normal-mode splitting function southwestern Weddell Sea, Geophys. J. Int., 203, 737–754. measurements, Geophys. J. Int., 184, 1223–1236. Wolstencroft, M., Shen, Z., Tor ¨ nqvist, T.E., Milne, G.A. & Kulp, M., 2014. Sasgen, I. et al., 2017. Joint inversion estimate of regional glacial isostatic Understanding subsidence in the Mississippi Delta region due to sediment, adjustment in Antarctica considering a lateral varying Earth structure ice, and ocean loading: Insights from geophysical modeling, J. geophys. (ESA STSE Project REGINA), Geophys. J. Int., 211, 1534–1553. Res., 119, 3838–3856. Schoof, C., 2007. Ice sheet grounding line dynamics: steady Wu, P., 1999. Modelling postglacial sea levels with power-law rheology and states, stability, and hysteresis, J. geophys. Res., 112, F03S28, a realistic ice model in the absence of ambient tectonic stress, Geophys. doi:10.1029/2006JF000664. J. Int., 139, 691–702. Schotman, H.H.A., Wu, P. & Vermeersen, L.L.A., 2008. Regional pertur- Wu, P., 2004. Using commercial ﬁnite element packages for the study of bations in a global background model of glacial isostasy, Phys. Earth earth deformations, sea levels and the state of stress, Geophys.J.Int., 158, planet. Inter., 171, 323–335. 401–408. Shapiro, N.M. & Ritzwoller, M.H., 2004. Inferring surface heat ﬂux dis- Wu, P. & Johnston, P., 1998. Validity of using ﬂat-earth ﬁnite element models tributions guided by a global seismic model: particular application to in the study of postglacial rebound, in Dynamics of the Ice Age Earth: a Antarctica, Earth planet. Sci. Lett., 223, 213–224. Modern Perspective, pp. 191–202, ed. Wu, P. Trans Tech Pub. Steffen, H., Kaufmann, G. & Lampe, R., 2014. Lithosphere and upper- Zhao, C., King, M.A., Watson, C.S., Barletta, V.R., Bordoni, A., Dell, M. mantle structure of the southern Baltic Sea estimated from modelling & Whitehouse, P.L., 2017. Rapid ice unloading in the Fleming Glacier relative sea-level data with glacial isostatic adjustment, Solid Earth, 5, region, southern Antarctic Peninsula, and its effect on bedrock uplift rates, 447–459. Earth planet. Sci. Lett., 473, 164–176. Steffen, H., Kaufmann, G. & Wu, P., 2006. Three-dimensional ﬁnite-element Zhao, S., Lambeck, K. & Lidberg, M., 2012. Lithosphere thickness and man- modeling of the glacial isostatic adjustment in Fennoscandia, Earth tle viscosity inverted from GPS-derived deformation rates in Fennoscan- planet. Sci. Lett., 250, 358–375. dia, Geophys. J. Int., 190, 278–292. van der Wal, W., Barnhoorn, A., Stocchi, P., Gradmann, S., Wu, P., Drury, Zhong, S., Paulson, A. & Wahr, J., 2003. Three-dimensional ﬁnite-element M. & Vermeersen, B., 2013. Glacial isostatic adjustment model with modelling of Earth’s viscoelastic deformation: effects of lateral variations composite 3-D Earth rheology for Fennoscandia, Geophys. J. Int., 194, in lithospheric thickness, Geophys. J. Int., 155, 679–695. 61–77. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press http://www.deepdyve.com/lp/oxford-university-press/the-impact-of-lateral-variations-in-lithospheric-thickness-on-glacial-nlNatm0257

Loading next page...

Publisher: Oxford University Press
Copyright: Copyright © 2022 The Royal Astronomical Society
ISSN: 0956-540X
eISSN: 1365-246X
DOI: 10.1093/gji/ggy158
Publisher site: See Article on Publisher Site

Abstract

Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Geophys. J. Int. (2018) 214, 811–824 doi: 10.1093/gji/ggy158 Advance Access publication 2018 April 19 GJI Geodynamics and tectonics The impact of lateral variations in lithospheric thickness on glacial isostatic adjustment in West Antarctica 1 1 2 2 Grace A. Nield, Pippa L. Whitehouse, Wouter van der Wal, Bas Blank, John 3 3 Paul O’Donnell and Graham W. Stuart Department of Geography, Durham University, Durham DH1 3LE, UK. E-mail: grace.a.nield@durham.ac.uk Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS, Delft, The Netherlands School of Earth and Environment, University of Leeds, Leeds LS2 9JT, UK Accepted 2018 April 19. Received 2018 April 17; in original form 2018 January 30 SUMMARY 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. Key words: Creep and deformation; Satellite geodesy; Antarctica; Dynamics of lithosphere and mantle; Rheology: crust and lithosphere; Rheology: mantle. due to large uncertainties that persist in both the ice-sheet history 1. INTRODUCTION since the Last Glacial Maximum (LGM) and the Earth structure in The process of Glacial Isostatic Adjustment (GIA) in Antarctica is this region. This has a direct impact on estimates of ice-mass loss well-studied (e.g. Whitehouse et al. 2012b;A et al. 2013; Argus derived from satellite gravimetry (e.g. the Gravity Recovery and et al. 2014) but GIA models continue to predict remarkably dif- Climate Experiment, GRACE) since Antarctic GIA is a signiﬁcant ferent present-day deformation rates (Mart´ ın-Espanol ˜ et al. 2016) The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 812 G. A. Nield et al. component of the total gravitational signal and must be removed to 2005; Nield et al. 2014; Wolstencroft et al. 2014). On the basis of yield estimates for ice-mass balance (King et al. 2012). wave speed variations, seismic studies can distinguish between ther- Traditionally, global and Antarctic-wide models of GIA have mal conduction and convection regimes in the upper mantle. The used a 1-D approximation of the Earth structure consisting of an conductive domain deﬁnes the tectonic plate. However, the elastic elastic lithosphere underlain by a linear viscoelastic upper and lower thickness varies as a function of timescale of surface loading and is mantle, where properties vary only radially (e.g. Peltier 1974; Milne typically thinner than the seismic lithospheric thickness. & Mitrovica 1996; Kendall et al. 2005). In reality, the structure of The studies and methods described above have used linear vis- the Earth is far more complex and models that reﬂect lateral as well coelastic rheology to model GIA. However, the use of power-law as vertical variations in Earth properties are needed to provide more rheology is becoming increasingly common (Wu 1999; Barnhoorn accurate predictions of present-day GIA-related deformation and et al. 2011; van der Wal et al. 2013; van der Wal et al. 2015). van geoid changes, both in Antarctica (A et al. 2013;van derWal et al. der Wal et al. (2015) used seismic velocity anomalies (Grand 2002) 2015;Sasgen et al. 2017) and elsewhere, for example, Greenland and geothermal heat ﬂux estimates (Shapiro & Ritzwoller 2004) (Khan et al. 2016). Including 3-D structure in GIA models is par- for Antarctica to infer mantle temperatures that were used to derive ticularly pertinent for Antarctica as this continent is considered to creep parameters for input to a power-law rheology. By deﬁning spa- consist of two distinct regions in terms of Earth structure: a thick tially varying creep parameters, the GIA model included laterally cratonic lithosphere and a high-viscosity uppermost mantle in the varying Earth structure. For this approach, the lithospheric thick- East, and thinner lithosphere and lower viscosity uppermost mantle ness is implicitly deﬁned by the creep parameters, rheological model in the West (Morelli & Danesi 2004). Modelling East and West and some threshold viscosity above which it can be considered to Antarctica with a 1-D Earth model as described above, therefore, behave elastically as described above. has the potential to produce incorrect estimates of the present- Modelling advances in the past few decades (Wu & Johnston day GIA signal in one or both of these regions. For example, A 1998; Latychev et al. 2005b;A et al. 2013;van derWal et al. 2013) et al. (2013) compared deformation rates predicted by a 3-D model have eased the computational burden of 3-D GIA modelling and incorporating laterally varying lithospheric thickness and mantle detailed data sets are now available that can be used to deﬁne lateral viscosity with a model that is the 1-D average of the 3-D proﬁle Earth structure (Ritsema et al. 2011;Heeszel et al. 2016). Hence, and found mismatches at Global Positioning System (GPS) loca- there are an increased number of studies incorporating 3-D Earth tions in Antarctica. Furthermore, capturing variability in the Earth structure into GIA models with both linear and nonlinear rheologies. structure within West Antarctica is important because regional 1-D Several approaches can be used to infer 3-D mantle viscosity (Ivins GIA studies have indicated differences in the Earth structure across & Sammis 1995;Kaufmann et al. 2005) and lithospheric thickness the region (Nield et al. 2014; Wolstencroft et al. 2015; Nield et al. for input to GIA models, with the latter being the focus of this 2016;Zhao et al. 2017). study. A seismically derived lithosphere–asthenosphere boundary This study focusses on how lateral variations in lithospheric thick- (LAB) depth is sometimes used to infer laterally varying GIA litho- ness impact predictions made by GIA models. The lithospheric spheric thickness with linear viscoelastic rheology, after scaling to thickness can be deﬁned by various criteria, such as a change in the account for differences between a seismically derived deﬁnition of method of heat transfer (Martinec & Wolf 2005), seismic anisotropy the lithosphere and the mechanical deﬁnition used in GIA studies. or resistivity (Eaton et al. 2009). In GIA modelling, the lithosphere For example, Kaufmann et al. (2000) reduced a seismically derived is deﬁned on the basis of mechanical properties and is considered LAB depth by a factor of two for their GIA modelling, and Khan to be a part of the crust and upper mantle that behaves elastically et al. (2016) used an adjustment factor to scale the LAB depth pub- on timescales of glacial cycles (tens of thousands of years; Mar- lished by Priestley & McKenzie (2013). However, it is not clear tinec & Wolf 2005; Watts et al. 2013; Kuchar & Milne 2015). The whether a lithosphere deﬁned by seismic properties can be con- lithosphere can be modelled with either a purely elastic rheology, verted to a lithosphere deﬁned by mechanical properties through a that is, it has no viscous component (e.g. Argus et al. 2014), or as a scaling factor. Seismic properties could have a different dependence viscoelastic material with sufﬁciently high viscosity that it does not on temperature and composition than mechanical properties. One way to circumvent this issue is to use temperatures derived from relax in response to surface loading on timescales of a glacial cycle (e.g. Kendall et al. 2005; Whitehouse et al. 2012b; Kuchar & Milne seismic velocity perturbations as input to a power-law rheology, 2015), thereby behaving elastically. Studies have also combined which eliminates the need to explicitly deﬁne lithospheric thick- these approaches, for example, Kaufmann et al. (2005) modelled a ness (van der Wal et al. 2013). In this approach, assumptions are 100 km thick lithosphere composed of a 30 km purely elastic layer made in converting seismic velocity anomalies into temperature and overlying a 70 km viscoelastic layer with a viscosity of 1 × 10 Pa viscosity, and the lithosphere is deﬁned implicitly by the effective s, which is approximately the limit of what could be considered viscosity. So although temperatures ultimately come from the same elastic over GIA timescales [e.g. 1 × 10 Pa s (Sasgen et al. 2017) seismic source as the LAB depths, no new assumptions are required –1 × 10 Pa s (Khan et al. 2016)]. Kuchar & Milne (2015)in- other than those made for converting seismic velocities to viscosity. vestigated the effect of depth-dependent viscosity in the lithosphere Previous studies investigating 3-D Earth structure in GIA mod- on relative sea-level predictions using a radially varying (i.e. 1-D) els of Antarctica have focussed on the effect of lateral variations Earth model and found that predictions made using a lithosphere in mantle viscosity (e.g. Kaufmann et al. 2005) or a combination with viscosity structure were similar to predictions made using a of laterally varying lithospheric thickness and upper-mantle viscos- purely elastic, but much thinner, lithosphere. ity (A et al. 2013; van der Wal et al. 2015) on present-day uplift To some extent, the apparent thickness of the lithosphere de- rates. Studies that isolate the effect of including lateral variations pends on the timescales of surface loading. Over long timescales in lithospheric thickness in models of GIA exist for regions in the (∼1 Myr), viscous relaxation in the lower lithosphere means that the northern hemisphere (Kaufmann et al. 2000; Zhong et al. 2003; lithosphere seems to behave as a relatively thin elastic layer (Watts Latychev et al. 2005a; Whitehouse et al. 2006;Steffen et al. 2014) et al. 2013). However, over GIA timescales (∼100 kyr), the litho- but currently no such study exists for Antarctica. sphere seems to behave as a thicker elastic layer (Martinec & Wolf Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 813 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 deﬁning a laterally varying lithospheric thickness mentioned above. The ﬁrst method (experiment 1) uses a scaled seismically inferred LAB depth to de- termine 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 differ- ences 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 litho- spheric 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 differ- ences 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 tem- perature as input to a power-law relationship, thereby assigning a viscosity structure to the lithosphere. In this method, the lithospheric thickness is implicitly deﬁned as the depth at which the resulting viscosity is small enough so that signiﬁcant deformation takes place Figure 1. Adjusted elastic lithospheric thickness derived from (a) Priestley during a glacial cycle. In order to determine the effect of includ- & McKenzie (2013)and (c)An et al. (2015a) LAB depths. Each colour ing viscosity in the lithosphere through a power-law rheology on represents a separate 20 km thick model layer with the lithospheric thickness being the upper bound of the colour, for example, orange denotes a LAB gradients of present-day uplift rate, we compare results using the depth/lithospheric thickness of 90 km. (b) and (d) show where the 3-D power-law rheology to those from the ﬁrst method which assumes lithosphere is thinner or thicker than the 90 km 1-D average. The regular that the lithosphere is elastic. For both methods, we use a ﬁnite mesh of 100 × 100 km is bounded by locations (−3000 km, −2000 km) element model with 20 km thick layers, meaning that the variable and (500 km, 2500 km), with an irregular mesh outside of this region. lithospheric thickness is captured in 20 km ‘steps’ between element locations. Due to large uncertainties in both Earth structure and at the edge of the model domain, for computational efﬁciency, and ice history, we do not attempt to ﬁt any observational data such as the domain has an overall width of 60 000 km. The extremely large GPS-observed uplift. model domain is required to ensure that boundary effects are negli- gible in the area of interest (Steffen et al. 2006). We do not model any ice loading changes outside Antarctica as the impact on uplift 2 METHODS AND DATA rates in Antarctica would be negligible (Whitehouse et al. 2012b) but we do include ocean loading. The model consists of 22 depth 2.1 Model and geometry layers representing the Earth’s surface to the core–mantle boundary (Table 1). A 30 km purely elastic crustal layer [the same thickness We use a 3-D ﬂat-earth ﬁnite element model constructed with the used by Kaufmann et al. (2005)] is underlain by eleven 20 km thick ABAQUS software package (Hibbitt et al. 2016) to compute the layers to 250 km depth to capture a higher resolution in the litho- solid Earth deformation in response to a changing surface load. sphere and uppermost mantle. Below this, layers are thicker (i.e. The validity of using this approach to model the Earth’s response lower resolution with depth) since surface deformation will be less to changes in ice-sheet loading has been shown previously by Wu sensitive to the details of mantle rheology below 250 km depth (Lau & Johnston (1998). This method has been used in many studies to et al. 2016). The buoyancy force is accounted for by applying Win- model GIA in regions such as Fennoscandia (Kaufmann et al. 2000; kler foundations to layer boundaries where a density contrast occurs Steffen et al. 2006), Antarctica (Kaufmann et al. 2005) and Iceland (Wu 2004). To ensure a direct comparability between the models, (Auriac et al. 2013), and has the advantage of computational efﬁ- the same mesh is used for all the experiments. ciency over a spherical global model. The ﬂat-earth ﬁnite element approach has been shown to be accurate when computing defor- mation within the ice margin for ice loads with comparable size to 2.2 Earth models and data the Laurentide Ice Sheet (Wu & Johnston 1998), which makes it applicable to the Antarctic Ice Sheet with its smaller lateral extent. The compressible elastic material properties for each layer described The mesh consists of eight-node brick elements with reduced above are listed in Table 1. The elastic and density structures of integration. The surface geometry of the mesh consists of a the Earth are derived from the preliminary reference Earth model 3500 × 4500 km area of interest embedded in a larger model (Dziewonski & Anderson 1981). For each element below the upper- domain. The area of interest represents West Antarctica and has most elastic layer, creep parameters are assigned on an element-by- elements of 100 × 100 km (elements are shown in Fig. 1). Outside element basis. The geometry of the mesh means that the lithospheric this region, the element size increases towards the periphery of the thickness (experiment 1) or the viscosity (experiment 2) can vary in model, from 550 km in East Antarctica to approximately 5000 km 20 km steps between the adjacent element locations. Latychev et al. Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 814 G. A. Nield et al. Table 1. Model layers and Earth parameters. Top of layer Top of layer Layer thickness Young’s modulus − 3 Layer radius (km) depth (km) (km) Density (kg m ) (GPa) Poisson’s ratio Rheology Lithosphere 6371 0 30 3196 173.9 0.28 Elastic Lithosphere or 6341 30 20 3379 173.9 0.28 Elastic/Power-law UM 6321 50 20 3377 173.9 0.28 lithosphere, or 6301 70 20 3375 173.3 0.28 linear viscoelastic 6281 90 20 3373 172.7 0.28 if UM 6261 110 20 3370 171.6 0.28 6241 130 20 3368 170.6 0.28 6221 150 20 3366 170.0 0.28 6201 170 20 3364 169.3 0.28 6181 190 20 3362 179.5 0.29 6161 210 20 3436 194.6 0.29 6141 230 20 3448 200.8 0.30 6121 250 80 3478 212.6 0.30 UM 6041 330 70 3525 224.4 0.30 Linear 5971 400 136 3812 277.2 0.29 viscoelastic- 5835 536 134 3978 377.8 0.28 variable LM 5701 670 251 4482 459.4 0.27 Linear viscoelastic 5450 921 250 4630 484.2 0.28 1 × 10 Pa s 5200 1171 430 4825 509.0 0.28 4770 1601 430 5036 570.1 0.29 4340 2031 430 5264 636.9 0.29 3910 2461 430 5464 704.5 0.30 (2005a) demonstrated that the differences in uplift rate over previ- Kaufmann et al. 2005; Kendall et al. 2005; Whitehouse et al. 2012b; ously glaciated regions in the northern hemisphere peak at 2 mm Ivins et al. 2013; Wolstencroft et al. 2015) and for the timescales −1 yr for a jump of 150 km between continental and oceanic litho- we are interested in the lithosphere is generally considered to be 24 25 spheric thickness, so we conclude that the effect of a 20 km jump elastic at viscosities above 1 × 10 –1 × 10 Pa s (Kaufmann et al. on the uplift rate is likely to be small. Our approach to deﬁning vari- 2005;Steffen et al. 2006; Barnhoorn et al. 2011;Khan et al. 2016). able lithospheric thickness allows us to use the same mesh for all Throughout the rest of this paper, we simply refer to this type of experiments, thus ensuring that the results are directly comparable. modelled lithosphere as the ‘elastic lithosphere’. For each model The sublithospheric upper mantle (down to a depth of 670 km) is a in experiment 1, we compare the laterally varying model with an linear viscoelastic layer with uniform viscosity and several different equivalent 1-D model in which the lithospheric thickness is simply upper-mantle viscosities are tested to determine dependence of the the average of the 3-D lithospheric thickness. results on the underlying viscosity (see Table2). Properties for the The second experiment uses a power-law rheology; this comple- lower mantle are the same for all models (Table 1). mentary approach allows us to investigate the differences in the The thickness of the lithosphere is deﬁned differently in ex- two methods used to deﬁne the lithospheric thickness in GIA mod- periments 1 and 2, as detailed in Table 2. The ﬁrst experiment elling. Mantle temperatures (An et al. 2015a) are used to determine considers an elastic lithosphere with spatial variability deﬁned by diffusion and dislocation creep parameters following the methods two different models of seismically derived LAB depth (Priestley described by Hirth & Kohlstedt (2003)and vander Wal et al. (2013; & McKenzie 2013;An et al. 2015a), as described in the follow- 2015). The reader is referred to these papers for a detailed descrip- ing sections. Seismically derived LAB depths tend to be thicker tion of the method. We limit the power-law rheology to the same than those inferred from GIA studies. For example, Steffen et al. horizontal and vertical domain as deﬁned by the An et al. (2015a) (2014) compare GIA-inferred lithospheric thicknesses in the Baltic LAB depths for two reasons: (1) so that the results of experiment 2 Sea with three LAB depth models and ﬁnd a consistently thinner can be directly compared with the results of experiment 1, and (2) lithosphere by 30–80 km. We, therefore, uniformly reduce the LAB so that the upper-mantle viscosity remains laterally uniform and, depths from the Priestley & McKenzie (2013)and An et al. (2015a) therefore, the effect of a spatially variable lithospheric thickness is models so that the thicknesses are more representative of a GIA isolated from all other parameters. 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 & 2.2.1 Priestley & McKenzie (2013) McKenzie (2013) because it is an Antarctic-speciﬁc model derived Priestley & McKenzie (2013) published a model of global seismic using many additional seismic stations. For those elements repre- velocities from the surface wave tomography which they used to senting the lithosphere, the viscosity is set to 1 × 10 Pa s to mimic derive mantle temperatures and lithospheric thickness on a 2 grid elastic behaviour on glacial timescales, apart from the uppermost with a resolution of 250 km horizontally and 50 km vertically. The 30 km layer which is modelled as purely elastic. This approach lithospheric thickness is deﬁned by the change in heat transfer from of modelling an elastic lithosphere using a viscoelastic rheology conduction to convection. We use this information in experiment 1a with high viscosity (including combinations of purely elastic and to deﬁne an elastic lithospheric thickness. The authors state that the viscoelastic rheology) is taken in many GIA studies (Peltier 2004; uncertainty on the lithospheric thickness is 20–30 km, and therefore Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 815 Table 2. Summary of the experiments and inputs used. Lithosphere Experiment deﬁnition Data used Ice models Upper-mantle viscosity (Pa s) Results 1a Elastic Priestley & McKenzie (2013) W12 (all viscosities); ICE-5G, 5 × 10 (all ice models); Comparison between 1-D and 19 20 21 LAB depths (adjusted) ICE-6G C, Uniform loading 5 × 10 ,1 × 10 ,1 × 10 3-D (for 5 × 10 Pa s) (W12 only) 1b Elastic An et al. (2015a) LAB depths W12 (all viscosities); ICE-5G, 5 × 10 (all ice models); Comparison between 1-D and 19 20 21 (adjusted) ICE-6G C, Uniform loading 5 × 10 ,1 × 10 ,1 × 10 3-D (for 5 × 10 Pa s) (W12 only) 2 Power-law Domain deﬁned by An et al. W12 5 × 10 Comparison between elastic (2015a) LAB depths, An et al. 3-D (1b) and power-law 3-D (2015b) temperatures used as (same ice model and input to power-law rheology upper-mantle viscosity) we reduce the LAB depths by 20 km to reﬂect the fact that a GIA- The LAB depths mapped onto the model elements is shown in inferred elastic lithosphere is typically thinner than a lithosphere Fig. 1(c), again with Fig. 1(d) showing where the 3-D model has based on the change in heat transfer method (Martinec & Wolf thicker or thinner lithosphere than the 1-D average [90 km, calcu- 2005; Priestley & McKenzie 2013). GIA-inferred LAB thicknesses lated over the Antarctic Plate which is the spatial limit of the inferred are less than the LAB depths in the Priestley & McKenzie (2013) LABdepthsinthe An et al. (2015a) model]. Similar to Priestley & model for regions such as Iceland [less than 50 per cent, 15–40 km McKenzie (2013), the lithosphere under East Antarctica is thicker (Auriac et al. 2013) compared with ∼95 km] and Fennoscandia than the 1-D average for the An et al. (2015a) model. The location [around 50–70 per cent, 93–110 km (Zhao et al. 2012) compared of the boundary between East and West Antarctica is similar in both with 120–200 km], and we use the uncertainty bound to reduce the seismic models, indicating that the uncertainty on the location of LAB depths to 70–90 per cent of the modelled values so that our the boundary is small. Some isolated regions of anomalously thick results are conservative. Fig. 1(a) shows the adjusted lithospheric lithosphere are also present in the Northern Antarctic Peninsula, thicknesses mapped onto the ABAQUS layers (Table 1), that is, at which the authors attribute to a remnant subducted slab from the each location on the mesh, the adjusted LAB depth is rounded to former subduction zone in this region. 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 2.3 Ice loading a thicker than average lithosphere. 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 2.2.2 An et al. (2015a,b) et al. 2013; Argus et al. 2014). Given this uncertainty, one of the The second model used in this study is from An et al. (2015a), who aims of this study is to investigate whether differences between 1- infer temperatures below Antarctica from the 3-D seismic velocity D and 3-D models are independent of the assumed ice history. To model AN1-S (An et al. 2015b), which has a horizontal resolution test this, we use several different deglacial models, or ice loading that increases from ∼120 km in the crust to ∼500 km at a depth of histories, that have quite different spatial patterns and magnitude 120 km and a vertical resolution of ∼25 to 150 km depth followed of loading changes, which, when applied to a speciﬁc Earth model, by ∼50 to 250 km depth. We use this model in two ways. First, the give different patterns of deformation. We compare gradients of seismically derived LAB, which is deﬁned by the depth where the present-day uplift rates between 1-D and 3-D Earth models using adiabat crosses the 1330 C geotherm, is used to deﬁne lithospheric the same ice history, revealing differences that may be directly at- thickness for the elastic lithosphere in experiment 1b. The uncer- tributed to the introduction of a varying lithospheric thickness, and tainty on the temperature is reported to be ±150 C, equivalent to then qualitatively compare the results from different ice models. ±15–30 km for the LAB depth, so we reduce the LAB depth by Three ice loading scenarios are used in the modelling: W12 15 km to be representative of GIA-elastic thicknesses, as explained (Whitehouse et al. 2012a), ICE-5G (Peltier 2004) and its successor in Section 1, and to be consistent with the scaling of the Priestley ICE-6G C (Argus et al. 2014; Peltier et al. 2015). The W12 ice & McKenzie (2013) model. Second, we use the temperatures di- loading model was derived using a glaciologically consistent nu- rectly as input to the power-law rheology in experiment 2. In this merical ice-sheet model that was tuned to provide the best possible experiment, we consider the 3-D spatial domain that deﬁnes the ﬁt to constraints of ice thickness change, whereas ICE-5G and ICE- lithosphere in experiment 1, but within this domain we use a power 6G C have been tuned to ﬁt ice thickness change constraints and law instead of elastic rheology. Results show the comparison of the GPS-observed uplift rates without satisfying ice-sheet physics. Fur- two 3-D models, thereby highlighting the differences due to rheo- thermore, in an attempt to fully isolate differences associated with logical deﬁnitions (see Table 2 for a summary of the models). In the introduction of a laterally varying lithospheric thickness from their model An et al. (2015a) do not infer temperatures for depths those caused by spatial variations in ice loading, we also construct shallower than 55 km. Therefore, when using the temperatures in an idealistic, spatially uniform loading history. In this scenario, the our model we specify a second elastic layer between 30 and 50 km amount of ice thickness change since the LGM is spatially uniform depth to bridge the gap between our uppermost elastic layer and the over the grounded area of the present-day ice sheet (as shown in temperature inputs. Table3 and Fig. 2) with a somewhat arbitrary 700 m of total ice loss Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 816 G. A. Nield et al. Table 3. Ice thickness change for the spatially uniform ice loading history for the West Antarctic Ice Sheet (WAIS) and the East Antarctic Ice Sheet (EAIS). Ice thickness change (m) Time period (ka) WAIS EAIS 20–15 −200 50 15–10 −300 60 10–5 −150 30 5–0 −50 10 Total: LGM to Present −700 150 Figure 3. Present-day uplift rates for the a) 1-D and b) 3-D models based on An et al. (2015a) LAB depths (experiment 1b), using the W12 ice loading history and upper-mantle viscosity 5 × 10 Pa s; c) difference in present- day uplift rates (1-D minus 3-D); d) difference in spatial gradient of uplift rate between 1-D and 3-D model (1-D minus 3-D) for the high resolution region of interest only—blue areas show where the 3-D model predicts higher amplitude and shorter wavelength deformation compared with the Figure 2. Regions of ice thickness change for the uniform loading his- 1-D model. tory for West Antarctica (blue) and East Antarctica (green). Key locations mentioned in the text also labelled. lithospheric thickness values to some of the models used in this study. However, we consider the impact of this to be small as there for West Antarctica and 150 m ice-sheet growth for East Antarc- −1 is, at most, ±0.7 mm yr difference to present-day uplift rates tica, applied over four time periods. These ice thickness changes when not including ocean loading at all; nevertheless, we choose approximate the mean ice loading changes in the W12 ice loading to include ocean loading with the intention of making the model as history [compare with ﬁg. 7 of Whitehouse et al. (2012a)]. realistic as possible. We keep the ocean loading the same for each ice model so that any differences in results may be attributed to dif- ferences in Earth properties. For the spatially uniform ice loading 2.4 Ocean loading history, we do not include any ocean loading since it is an idealized The model approach we have used in this study does not solve loading history and would not produce a realistic sea-level change the sea-level equation (Farrell & Clark 1976) and cannot compute when modelled with a global GIA model. variable sea level with time. We, therefore, take the approach of applying an ocean load that has been derived using a global, spher- ically symmetric GIA model (Mitrovica & Milne 2003; Kendall 3 RESULTS et al. 2005; Mitrovica et al. 2005). The GIA model uses a given ice loading history and an Earth model to calculate changes in In order to determine the effect of introducing a varying lithospheric sea level (i.e. a change in surface loading due to a change in the thickness in experiment 1 (elastic case), we examine the differences depth of the ocean) with time. The ocean load was computed using between the 1-D and 3-D model output in terms of the spatial the ice loading histories W12, ICE-5G and ICE-6G C in combina- gradient of predicted present-day uplift rates. The spatial gradient tion with a three-layer Earth model and the output is a time- and is simply the derivative of the present-day uplift rate ﬁeld and we space-variable load that can be applied to our laterally varying ﬂat take the scalar magnitude of the gradient (i.e. it is always positive) Earth model. We use an Earth structure that is representative of since the direction of the slope is not of interest. We calculate the our 1-D average models with a lithospheric thickness of 96 km, spatial gradient over the 100 km resolution area of interest only. upper-mantle viscosity of 5 × 10 Pa s, and lower-mantle viscosity Differences in present-day uplift rates are relatively small 22 −1 of 1 × 10 Pa s. We acknowledge the inconsistencies inherent in (± 3mmyr , Fig. 3c) and the sign of the difference does not this approach in that the ocean load is computed using a 1-D Earth yield useful information. For example, in the Siple Coast the 3-D model that may have different average upper-mantle viscosity and model predicts greater uplift at the coast but also more subsidence Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 817 in the interior of West Antarctica, in other words, the 1-D model ICE-6G C ice model (Figs 5c and g) shows a prominent blue bull’s- under-predicts the magnitude of the response compared with the eye located near the Siple Coast related to a large unloading event. 3-D model (Figs 3a and b). Differencing the deformation rates (1-D The unloading results in steeper uplift gradients and a higher peak minus 3-D) shows both a positive and negative difference (Fig. 3c), amplitude in the 3-D case compared with the 1-D case, but in the masking the fact that the 3-D model produces a higher peak-to-peak centre, that is, at the peak itself, the gradient is the same (white difference in uplift rate (i.e. higher peaks and lower troughs). Calcu- on the ﬁgures). This effect can also be observed with the spatially lating the spatial gradient of uplift rate for the 1-D and 3-D models uniform loading history (Figs 5d and h), where the periphery of and differencing them indicates how the amplitude and wavelength the ice sheet shows the most sensitivity to variations in lithospheric of deformation varies between the two models (Fig. 3d). A higher thickness (i.e. largest differences in predicted present-day uplift rate amplitude and shorter wavelength response would be expected from gradients), and the interior shows little difference between the 1-D a thinner lithosphere compared with a thicker lithosphere (Wolsten- averaged model and the 3-D model. croft et al. 2015). This canbeobservedinFig. 4 from the proﬁle of Despite the localized differences in spatial pattern, all combi- uplift rate and gradient of uplift rate across the Antarctic Peninsula, nations of ice loading history and LAB model tested here yield where the lithospheric thickness in the 3-D model is thinner than the same ﬁrst-order result across most of West Antarctica—use of that of the 1-D average. The uplift rate predicted by the 3-D model a 1-D averaged lithospheric thickness results in lower magnitude (orange solid line in Fig. 4) has a higher amplitude and shorter gradients (lower amplitude and longer wavelength) of present-day wavelength (by one grid cell) than the 1-D model (green solid line uplift rate compared with the 3-D case, and hence predominantly in Fig. 4). This means that the gradient of uplift in the 3-D model negative differences across West Antarctica in Fig. 5. Any positive will be steeper around the peak of the rebound, as indicated by the (red) differences in West Antarctica result from the longer wave- blue colour in the inset, but it tails off more quickly than in the 1-D length deformation predicted by the 1-D model resulting in steeper model resulting in the 1-D gradient being steeper at the periphery gradients than the 3-D model at the periphery of the rebound. This (indicated by red in the inset). This gives a characteristic pattern result is insensitive to the ice model used, although the actual spa- of a white bulls eye (where the gradients are the same at the tip tial patterns shown in Fig. 5 do depend on the ice loading history of the peak), surrounded by blue where the 3-D gradient is higher since the biggest differences in gradients when comparing a uniform (negative gradient difference), with red at the periphery (Fig. 3d). lithospheric thickness to a laterally varying lithospheric thickness In East Antarctica, where the 1-D averaged lithospheric thickness mostly occur around the margins of loading/unloading centres. The is thinner than the 3-D model, and the present-day uplift rate gra- ice loading histories used in this study neglect any changes in ice- dients are steeper in the 1-D model output, the gradient difference sheet thickness over the past few thousand years, such as those is positive and shown as red, with the same characteristic white observed in the Antarctic Peninsula (Nield et al. 2012) and Siple at the peak of the uplift/subsidence centres (Fig. 3d). Results for Coast (Catania et al. 2012). Including these late Holocene changes experiment 1 in Sections 3.1–3.3 are shown in the same format as would have the effect of changing the pattern of localized differ- Fig. 3(d)—as differences in uplift rate gradient between 1-D and ences providing the underlying mantle viscosity was sufﬁciently 3-D models for our high resolution region of interest. low to respond on a timescale of ∼2000 yr or less. 3.2 Effect of LAB model 3.1 Effect of ice loading history The choice of LAB model used to deﬁne spatial variations in litho- Fig. 5 shows the difference in spatial gradient of the present-day spheric thickness has the potential to inﬂuence the results. The An uplift rate when comparing the 3-D and 1-D models for the four et al. (2015a) model has a higher resolution than the Priestley & different ice loading histories used in this study. Results are shown McKenzie (2013) model and therefore contains more spatial vari- for both models of lithospheric thickness used in experiment 1 – ability in the LAB depths. The bottom row of Figs 5 and 6 show Priestley & McKenzie (2013)and An et al. (2015a); the bottom the difference in uplift rate between the two LAB models, for the row in Fig. 5 shows the difference in uplift rate between these two different ice loading models and upper-mantle viscosities, respec- models. The upper- and lower-mantle viscosities are kept the same tively. The impact of the LAB model in isolation can most clearly 20 22 for all models (5 × 10 and 1 × 10 Pa s, respectively). This plot be observed in Fig. 5l—the model that uses the uniform loading can help us to understand what effect the ice loading history has on history—because there are no spatial variations in ice loading that the results. can amplify signals. The differences in Fig. 5l directly reﬂect the Each ice loading history results in different localized spatial pat- differences between the two LAB models (Fig. 1), with the greatest terns of present-day uplift rate gradient reﬂecting the spatial vari- effects being in the Northern Antarctic Peninsula, where An et al. ability of ice loading or unloading between the models. Differences (2015a) identify a region of anomalously thick lithosphere, and in in the present-day uplift rate gradients of the 3-D and 1-D mod- Coats Land (Fig. 2) where the boundary between East and West els, whether negative or positive, are focussed around the margins Antarctica is deﬁned differently for each model. The differences −1 of centres (or ‘peaks’) of uplift or subsidence. This is because peak at ±3.5 mm yr for this latter region when using the W12 ice the lithospheric thickness in the 3-D model is thinner/thicker than loading history (Fig. 5i) because large loading changes across this the 1-D average and hence produces higher/lower amplitude and region during the past 5 ka (Whitehouse et al. 2012a) amplify the steeper/shallower gradients than the 1-D model. When comparing signal. All other ice loading/mantle viscosity combinations result −1 peaks of uplift rate gradient between the 1-D and 3-D model (for in differences of ±1.5 mm yr or less. the same ice history) they have different amplitudes but the gradi- We can draw several conclusions from these observations. First, ents at the crest of the peaks will often be the same or very similar, the results are more dependent on the ice loading history used than as explained previously, resulting in a small area of white at the the choice of LAB model. Second, we do not gain signiﬁcant extra centre of the region of uplift/subsidence (Fig. 5). For example, the information by using a higher resolution LAB model that resolves Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 818 G. A. Nield et al. Figure 4. Uplift rate (left-hand axis) for the 1-D (solid green) and 3-D (solid orange) models along the proﬁle shown in the inset. Also shown is the gradient of uplift rate (right hand axis) along the proﬁle for the 1-D (dashed green) and 3-D (dashed orange) models, with shading according to the difference in gradient shown in the inset (1-D minus 3-D; same as Fig. 3d). Black dashed line indicates the difference in gradient shown in the inset plot. Figure 5. Difference in spatial gradient of uplift rate between 1-D and 3-D models (1-D minus 3-D) for ice loading histories (from left to right) W12 (a, e), ICE-5G (b, f), ICE-6G C (c, g) and the uniform loading history (d, h), and for the two different LAB models, Priestley & McKenzie (2013;top row) andAn et al. (2015a; middle row). All models have an upper-mantle viscosity of 5 × 10 Pa s. The dashed-dotted black line delineates where the 3-D lithosphere is thinner or thicker than in the 1-D case, as shown in Figs 1(b) and (d). Panels (i)-(l) show the difference in uplift rate between the 3-D LAB models (Priestley & McKenzie (2013) minus An et al. (2015a)). smaller scale variations in lithospheric thickness, even if we increase both seismically inferred LAB models show a clear East–West di- the GIA model horizontal resolution to 50 km (Section 4.4). Finally, vide, with the East having thicker-than-1-D-average lithosphere and Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 819 Figure 6. Difference in spatial gradient of uplift rate between 1-D and 3-D models (1-D minus 3-D) for different values of upper-mantle viscosity (from left to right), for the two different LAB models, Priestley & McKenzie (2013;top row) andAn et al. (2015a; middle row), and using the W12 ice history. The dashed-dotted black line delineates where the 3-D lithosphere is thinner or thicker than in the 1-D case, as shown in Figs 1(b) and (d). Panels (i)–(l) show the difference in uplift rate between the 3-D models for the two different LAB models [Priestley & McKenzie (2013) minus An et al. (2015a)]. the West having thinner-than-1-D-average lithosphere, as indicated the western Ross Sea between 10 and 5 ka, whereas rebound in the by the dashed-dotted lines in Figs 5 and 6. This demarcation co- lower viscosity models (Figs 6a and b, e and f) is dominated by incides with the regions where the amplitude of gradients of uplift the response to late Holocene ice thinning along the Siple Coast rates for the 3-D model are higher (in West Antarctica) or lower (in and Southern Antarctic Peninsula, as indicated by the blue areas in East Antarctica) than the 1-D model and it is clearly the feature that Figs 6(a and b) and (e and f). has the most impact on gradients of uplift rates. 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 viscosi- ties we observe the same result of higher amplitude and shorter 3.3 Effect of upper-mantle viscosity wavelength deformation in West Antarctica for the 3-D model (blue Upper-mantle viscosity exerts a strong control on mantle relaxation in the ﬁgures) supporting the hypothesis that the lithospheric thick- times and hence uplift rates. To test whether our results are de- ness controls the wavelength of the signal captured in the modelling. pendent on the underlying upper-mantle viscosity, we calculated the difference in present-day uplift rate gradients using four upper- mantle viscosities, for both the LAB models in experiment 1, using 3.4 Effect of power-law rheology in the lithosphere just the W12 ice loading history (Fig. 6). Comparing the results, we see similar patterns of gradient dif- Modelling the lithosphere using a power-law rheology means that ferences for the weaker upper-mantle viscosities (5 × 10 and there is the potential for it to deform viscously, depending on the 1 × 10 Pa s, Figs 6a and b, e and f) and the stronger upper-mantle input temperature used to derive creep parameters and the stress 20 21 viscosities (5 × 10 and 1 × 10 , Figs 6c and d, g and h), although from the ice loading. We compare results using power-law rheology the two sets of patterns are different from each other. The two sets of (experiment 2) and input temperatures from An et al. (2015a; Fig. 7) patterns reﬂect sensitivity to different periods in the deglacial his- with results from the equivalent experiment 1b model that has a spa- tory of the W12 ice model (Whitehouse et al. 2012a). The models tially variable elastic lithosphere (Fig. 8); the two models have the with stronger mantle viscosities and slower relaxation time (Figs 6c same laterally varying lithospheric thickness but different rheology and d, g and h) are still rebounding in response to ice thinning in (see also Table 2). The upper-mantle viscosity (5 × 10 Pa s) and Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 820 G. A. Nield et al. Figure 7. Top row: temperatures from the An et al. (2015b) model averaged over the ﬁnite element model layers. Bottom row: effective viscosity at the present-day for the same model layers as the top row, calculated following the methods detailed in van der Wal et al. (2015). Red circle in panel (d) shows low viscosity lithosphere mentioned in the text. Elements below the spatially variable lithospheric thickness from An et al. (2015a) are greyed out (cf. Fig. 1c). Figure 8. (a) Difference in spatial gradient between the 3-D elastic-only case (experiment 1b) and the 3-D power-law case (experiment 2; elastic-only case minus power-law case), for the W12 ice loading history with upper-mantle viscosity of 5 × 10 Pa s. (b) Proﬁle of uplift rate for the elastic (green solid) and power-law (orange solid) cases and the gradient of each (dashed lines, right hand axis) along the proﬁle shown in (a). ice loading history (W12) are the same for both models. Modelling The patterns of gradient difference show in Fig. 8(a) are unlike the lithosphere with a power-law rheology has the effect of reducing the previous results. Around the Weddell Sea (Fig. 2), there is a dark the local effective elastic thickness (cf. Kuchar & Milne 2015)sowe blue region indicating higher amplitude deformation in the power- expect the power-law lithosphere (experiment 2) to behave as if it law model compared with the elastic model in experiment 1, which were thinner than the elastic lithosphere (experiment 1). In Fig. 8(a), may be related to the relatively low viscosity in the lithosphere at we plot the difference in uplift rate gradient as elastic minus power- 70–90 km depth (around 1 × 10 Pa s, see the red circle in Fig. 7d) law so that the colour scale can be compared with the earlier plots compared with the elastic lithosphere case (1 × 10 Pa s). Since the of 1-D minus 3-D. The effective viscosities for elements that lie in viscosity is dependent on the stress induced from ice load changes, the lithosphere are also calculated, following the methods described the low viscosity in this region may also be associated with late ice in van der Wal et al. (2015), and shown in Fig. 7 along with the loading changes deﬁned within the W12 model. In fact, viscosity temperatures from the An et al. (2015a) model that were used to in this region is up to an order of magnitude lower (1 × 10 Pa derive the creep parameters. 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/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 821 1-D versus 3-D is no longer present. This may be related to the because material properties determine the viscosity depending on fact that in this region, the seismic data indicate that there is a slab timescale and this would, for example, allow relaxation of the lower of relatively cold material at a depth of 50–70 km, resulting in a lithosphere over multiple glacial cycles. relatively high viscosity and therefore a very similar response to the It is, therefore, important to consider both how the lithosphere case with the elastic lithosphere in experiment 1. is deﬁned and how thickness variations are accounted for in the The proﬁles of present-day uplift rate and uplift rate gradient next generation of 3-D GIA models. As a minimum, East and West shown in Fig. 8(b) demonstrate that in the experiment that uses Antarctica should be considered separately in terms of Earth struc- power-law rheology, the peaks have a higher amplitude and shorter ture as both seismically derived LAB models considered here show wavelength than in the elastic lithosphere experiment. For the 50– a clear East-West divide in lithospheric thickness. We have shown 70 km depth layer, the viscosity within the lower lithosphere under that a model with higher resolution spatial variability in lithospheric 20–21 West Antarctica is around 1 × 10 Pa s, meaning it will de- thickness makes little difference to our results, however, represent- form viscously on glacial timescales of tens of thousands of years. ing lithospheric thickness variations within West Antarctica will This means that when using a power-law rheology to model the become more important as ice loading histories evolve to contain lithosphere only the upper 50 km will behave elastically over the greater spatial detail and include changes in ice thickness over the timescales of interest (cf. Kuchar & Milne 2015). 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 predic- tions of uplift rate and ultimately a more accurate GIA correction for 4 DISCUSSION GRACE data. This is particularly pertinent for the dynamic region of West Antarctica that is currently experiencing a large amount of 4.1 Implications for future GIA models ice-mass loss (Rignot et al. 2014). 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 4.2 Implications for interpretation of observations of GIA in higher gradients of predicted present-day uplift rates (i.e. higher Geodetic observations of bedrock deformation provide useful data amplitude and shorter wavelength) in West Antarctica compared with which to constrain models of GIA. Consideration of laterally with a uniform elastic lithospheric thickness that is simply the av- varying Earth structure may result in a better ﬁt between model erage of the former. We have made this comparison, ﬁrst of all, to predictions and observations in some areas. For example, Wolsten- isolate the effect of introducing variable lithospheric thickness from croft et al. (2015) could not ﬁt the spatial pattern of GPS-observed any other factors that perturb predictions of uplift rates, and second, uplift in the southern Antarctic Peninsula with a 1-D Earth struc- because many global GIA models use a 1-D Earth model derived ture having tested several variations on recent ice loading history. from globally averaged parameters. The mean lithospheric thickness It is possible that the strong spatial gradient in uplift revealed by over the GIA model domain of both models of seismically derived differencing GPS rates recorded at sites on the east and west of the LAB depth used in experiment 1 is 90 km, similar to values used Antarctic Peninsula could be explained with the introduction of a in studies of global GIA [80–90 km, Mitrovica & Forte (2004)and thinner lithosphere in this region (e.g. 50–70 km, Fig. 1), which Peltier et al. (2015), respectively]. Our results indicate that global would be able to capture shorter wavelength differences in uplift, as 1-D GIA models with a ∼90 km lithospheric thickness would pre- we have shown. However, before such a comparison is made, Late dict lower amplitude and longer wavelength uplift rates across West Holocene ice loading changes (e.g. Nield et al. 2012; Nield et al. Antarctica than would be predicted with a more realistic, spatially 2016) must be incorporated into current deglacial models. variable lithosphere. This means that uplift rates, and hence geoid Future observations of GIA should aim to be positioned in loca- changes, would be smoothed out over a wider area potentially lead- tions that would help to constrain 3-D Earth structure. In particular, ing to an inaccurate GIA correction for GRACE data. A 1-D model increasing the density of GPS networks across West Antarctica with a lithospheric thickness representative of the average of West would provide additional constraints for determining lithospheric Antarctica (70 km) produces a closer match to results from the 3-D thickness because shorter wavelength solid Earth deformation could model than the 1-D Antarctic average lithosphere (90 km), apart be observed. For example, Nield et al. (2014)wereabletomore from in regions where the lithosphere is even thinner (e.g. Southern tightly constrain lithospheric thickness in the northern Antarctic Antarctic Peninsula, 50 km, Fig. 1c). This suggests that modelling Peninsula by using observations from the dense LARISSA network. East and West Antarctica with a separate 1-D Earth model is an Furthermore, measurements along the boundary between East and important ﬁrst step in improving GIA models of Antarctica. West Antarctica would provide useful information in delimiting Furthermore, modelling the lithosphere with a power-law rhe- this transition in Earth structure for the purposes of GIA models. ology has the effect of reducing the thickness of the GIA litho- Additional measurements of horizontal deformation could be in- sphere (i.e. the portion of the lithosphere acting elastically on GIA strumental in constraining lateral variations in Earth structure in timescales) compared with the elastic case because the viscosity this region. 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 4.3 Implications for ice-sheet models experiment 2, we have shown that using these two different deﬁni- tions of the lithosphere leads to differences in gradients of present- We have demonstrated that the areas most affected by the inclusion day uplift rates despite input parameters (i.e. seismically derived of a spatially variable lithospheric thickness lie around the margins LAB depth and seismically derived temperatures) ultimately com- of ice loading changes, including (for most combinations of ice his- ing from the same source. Using a power-law rheology provides tory and Earth model tested) the Amundsen Sea sector and the Siple a more consistent way of modelling GIA over multiple timescales Coast (locations shown in Fig. 2). This has important implications Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 822 G. A. Nield et al. for ice dynamics in marine-grounded areas that lie on a reverse thickness when modelling the solid Earth response to surface load- slope bed, for example, West Antarctica. Grounding line dynamics ing across West Antarctica. The ﬁrst experiment used estimates for control ice-sheet stability and evidence shows that a reverse slope the depth of the LAB derived from seismic studies to model the bed can reduce ice-sheet stability because, as the grounding line lithosphere as an elastic layer, an approach taken in many GIA retreats into deeper water, ice ﬂux across the grounding line will studies. We have compared results from 3-D models (varying litho- increase, potentially leading to net ice loss and hence further retreat spheric thickness) and equivalent 1-D models (uniform lithospheric (e.g. Schoof 2007). thickness that is the average of the 3-D model). For all combinations Studies of Antarctic ice loss that make use of a coupled ice- of ice history, LAB model and underlying upper-mantle viscosity sheet–sea-level model have shown that bedrock uplift has a stabi- tested, we ﬁnd that the use of a 1-D averaged lithospheric thickness lizing effect on marine-grounded ice due to reducing the slope of a results in lower gradients (i.e. lower amplitude and longer wave- reverse bed, resulting in less ice loss from Antarctica (Gomez et al. length) of uplift rate compared with the use of a spatially variable 2010; Gomez et al. 2013). Including a spatially variable lithospheric (thinner in West Antarctica) lithospheric thickness. This means that thickness would increase the stabilizing effect of bedrock uplift on the present-day uplift rate is smoothed over a wider area in the 1-D the marine-grounded sector of the ice sheet in West Antarctica com- model and the magnitude of peaks and troughs of deformation is pared with a 1-D averaged model because, as we have shown, the smaller. This has important implications for ice-sheet modelling thinner lithosphere results in higher amplitude uplift in the inte- studies as steeper spatial gradients of uplift may promote stability rior, thereby reducing the slope of the reverse bed further. This in marine-grounded regions of West Antarctica. has been demonstrated by Gomez et al. (2015) and Pollard et al. The biggest difference in results between the two different seis- (2017) who show that a 1-D Earth model with a 50 km lithospheric mically derived LAB models used is in the Northern Antarctic thickness and low mantle viscosity results in increased stabiliza- Peninsula and at the boundary between East and West Antarctica, tion over a 1-D model with thicker lithosphere and higher mantle partly due to the An et al. (2015a) model having higher resolution viscosity. Furthermore, Gomez et al. (2018) showed that a cou- and a greater level of detail. The most important feature of these pled ice-sheet–sea-level model with a 3-D Earth structure (laterally LAB models is the delineation of where the lithosphere is thinner varying lithospheric thickness and upper-mantle viscosity) results than average in West Antarctica, which is a stable feature across in signiﬁcant regional differences in ice-sheet thickness when com- different seismic models, although the location of this boundary pared with results using a 1-D Earth structure. In particular, their is important as it can affect uplift rates in this area. Within West model predicts thicker ice and less retreat of the grounding line Antarctica the localized patterns of differences in uplift rate gradi- over the last deglaciation at the periphery of the Ross Sea region ent are sensitive to the choice of ice loading history, with largest (Fig. 2) where the lithosphere is thinner, and upper-mantle viscosity differences focussed around centres of loading or unloading. The is lower, than their 1-D average model. Including 3-D Earth struc- choice of underlying mantle viscosity also plays a role because ture in GIA models and ice dynamic models is, therefore, necessary the viscosity deﬁnes the relaxation time of the mantle, which in for determining the dynamics of past ice-sheet change and accu- turn determines which regions will still be deforming in response rately assessing the current and future state of the West Antarctic to past ice-sheet change. Including a laterally variable lithospheric Ice Sheet. thickness within West Antarctica will become even more important once ice loading histories incorporate changes from the past few thousand years. 4.4 Limitations The second experiment in this study investigated the difference between two methods of deﬁning the lithosphere in GIA modelling. Model resolution is an important consideration for any GIA model. We compared the elastic lithosphere in experiment 1 with the use Here, we restricted the spatial resolution to 100 × 100 km elements of power-law rheology in experiment 2, which deﬁnes viscosity in the area of interest, purely for computational efﬁciency. We tested based on material parameters and loading changes, and hence im- the effect of running a higher resolution model, increasing the mesh plicitly deﬁnes the lithosphere based on whether the viscosity is resolution to 50 km in the area of interest. While the output is high enough to behave elastically over the timescale in question. smoother, the 50 km resolution model did not reveal any additional Our results demonstrate that using a power-law rheology produces features that are not captured by the 100 km mesh and considering higher amplitude peaks of deformation than using a 3-D elastic- the extra computation time, we conclude that the coarser resolution only lithosphere because in the power-law case, the thickness of the is satisfactory for the experiments performed in this study. portion of the lithosphere that behaves elastically is reduced. Deﬁn- In the ﬁnite element model, material properties are considered ing the lithosphere in this way could provide a more robust model compressible in the computation of deformation, but the effect this of GIA since the thickness of the lithosphere is less rigidly deﬁned has on buoyancy forces is not considered. The model also neglects than in the elastic (i.e. very high viscosity) case and relaxation in self-gravitation, that is, changes in gravitational potential caused by the lower lithosphere could be important when modelling several deformation, which is a feature of most spherical models. However, glacial cycles (Kaufmann et al. 2005). Schotman et al. (2008) state that when using a ﬂat-earth model, Future GIA models should seek to include a spatially vary- the lack of sphericity partly cancels the lack of self-gravitation. ing lithospheric thickness, or at the very least to represent thin- Furthermore, since we are looking at differences between models, ner/thicker lithosphere in West/East Antarctica; we ﬁnd that inclu- any errors arising due to the neglect of such features will effectively sion of this transition has a ﬁrst order effect on the predicted pattern be cancelled out. of present-day deformation. Regional 1-D GIA models should en- sure that the local lithospheric thickness is adequately represented rather than using an average of a wider Antarctic domain. Fur- 5 CONCLUSIONS thermore, including a spatially variable lithosphere could lead to a We have presented the results of two experiments that seek to in- better ﬁt to GPS-observed uplift rates, especially in regions where a vestigate the impact of including lateral variations in lithospheric thinner lithosphere might be necessary to capture shore wavelength Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 Impact of laterally varying lithosphere 823 signals. This, in turn, could improve GIA models in West Antarctica Heeszel, D.S. et al., 2016. Upper mantle structure of central and West Antarc- tica from array analysis of Rayleigh wave phase velocities, J. geophys. where the uncertainty is large, although lateral variations in mantle Res., 121, 1758–1775. viscosity and better constraints on ice history would also be required Hibbitt, D., Karlsson, B. & Sorensen, P., 2016. Getting Started with ABAQUS, to provide an improved correction for GRACE data. Version (6.14), Hibbitt, Karlsson & Sorensen, Inc. Hirth, G. & Kohlstedt, D., 2003. Rheology of the upper mantle and the mantle wedge: a view from the experimentalists, in Inside the Subduction ACKNOWLEDGEMENTS Factory, ed. Eiler, J., pp. 83–105, American Geophysical Union. Ivins, E.R., James, T.S., Wahr, J., Schrama, E.J.O., Landerer, F.W. & Simon, We thank Rebekka Steffen and an anonymous reviewer for their K.M., 2013. Antarctic contribution to sea level rise observed by GRACE constructive comments that helped to improve the manuscript. with improved GIA correction, J. geophys. Res., 118, 3126–3141. GAN and JPOD are supported by NERC grant NE/L006294/1. Ivins, E.R. & Sammis, C.G., 1995. On lateral viscosity contrast in the mantle PLW is a recipient of an NERC Independent Research Fellow- and the rheology of low-frequency geodynamics, Geophys. J. Int., 123, ship (NE/K009958/1). This research is a contribution to the SCAR 305–322. SERCE program. All ﬁgures have been produced using the GMT Kaufmann, G., Wu, P. & Ivins, E.R., 2005. Lateral viscosity variations software package (Wessel & Smith 1998). beneath Antarctica and their implications on regional rebound motions and seismotectonics, J. Geodyn., 39, 165–181. Kaufmann, G., Wu, P. & Li, G., 2000. Glacial isostatic adjustment in Fennoscandia for a laterally heterogeneous earth, Geophys. J. Int., 143, REFERENCES 262–273. A, G, Wahr, J. & Zhong, S., 2013. Computations of the viscoelastic response Kendall, R.A., Mitrovica, J.X. & Milne, G.A., 2005. On post-glacial sea of a 3-D compressible Earth to surface loading: an application to Glacial level – II. Numerical formulation and comparative results on spherically Isostatic Adjustment in Antarctica and Canada, Geophys. J. Int., 192, symmetric models, Geophys. J. Int., 161, 679–706. 557–572. Khan, S.A. et al., 2016. Geodetic measurements reveal similarities between An, M. et al., 2015a. Temperature, lithosphere–asthenosphere boundary, and post–last glacial maximum and present-day mass loss from the Greenland heat ﬂux beneath the Antarctic Plate inferred from seismic velocities, J. ice sheet, Sci. Adv., 2, doi:10.1126/sciadv.1600931. geophys. Res., 120, 8720–8742. King, M.A., Bingham, R.J., Moore, P., Whitehouse, P.L., Bentley, M.J. & An, M. et al., 2015b. S-velocity model and inferred Moho topography Milne, G.A., 2012. Lower satellite-gravimetry estimates of Antarctic sea- beneath the antarctic plate from rayleigh waves, J. geophys. Res., 120, level contribution, Nature, 491, 586–589. 359–383. Kuchar, J. & Milne, G.A., 2015. The inﬂuence of viscosity structure in the Argus, D.F., Peltier, W.R., Drummond, R. & Moore, A.W., 2014. The Antarc- lithosphere on predictions from models of glacial isostatic adjustment, J. tica component of postglacial rebound model ICE-6G C (VM5a) based Geodyn., 86, 1–9. on GPS positioning, exposure age dating of ice thicknesses, and relative Latychev, K., Mitrovica, J.X., Tamisiea, M.E., Tromp, J. & Moucha, R., sea level histories, Geophys. J. Int., 198, 537–563 2005a. Inﬂuence of lithospheric thickness variations on 3-D crustal veloc- Auriac, A., Spaans, K.H., Sigmundsson, F., Hooper, A., Schmidt, P. & Lund, ities due to glacial isostatic adjustment, Geophys. Res. Lett., 32, L01304. B., 2013. Iceland rising: solid Earth response to ice retreat inferred from Latychev, K., Mitrovica, J.X., Tromp, J., Tamisiea, M.E., Komatitsch, D. & satellite radar interferometry and visocelastic modeling, J. geophys. Res., Christara, C.C., 2005b. Glacial isostatic adjustment on 3-D Earth models: 118, 1331–1344. a ﬁnite-volume formulation, Geophys. J. Int., 161, 421–444. Barnhoorn, A., van der Wal, W., Vermeersen, B.L.A. & Drury, M.R., 2011. Lau, H.C.P., Mitrovica, J.X., Austermann, J., Crawford, O., Al-Attar, D. & Lateral, radial, and temporal variations in upper mantle viscosity and Latychev, K., 2016. Inferences of mantle viscosity based on ice age data rheology under Scandinavia, Geochem. Geophys. Geosyst., 12, Q01007, sets: Radial structure, J. geophys. Res., 121, 6991–7012. doi:10.1029/2010GC003290. Mart´ ın-Espanol, ˜ A., King, M.A., Zammit-Mangion, A., Andrews, S.B., Briggs, R.D. & Tarasov, L., 2013. How to evaluate model-derived deglacia- Moore, P. & Bamber, J.L., 2016. An assessment of forward and inverse tion chronologies: a case study using Antarctica, Quat. Sci. Rev., 63, GIA solutions for Antarctica, J. geophys. Res., 121, 6947–6965. 109–127. Martinec, Z. & Wolf, D., 2005. Inverting the Fennoscandian relaxation- Catania, G., Hulbe, C., Conway, H., Scambos, T.A. & Raymond, C.F., 2012. time spectrum in terms of an axisymmetric viscosity distribution with a Variability in the mass ﬂux of the Ross ice streams, West Antarctica, over lithospheric root, J. Geodyn., 39, 143–163. the last millennium, J. Glaciol., 58, 741–752. Milne, G.A. & Mitrovica, J.X., 1996. Postglacial sea-level change on a Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference Earth rotating Earth: ﬁrst results from a gravitationally self-consistent sea-level model, Phys. Earth planet. Inter., 25, 297–356. equation, Geophys. J. Int., 126, F13–F20. Eaton, D.W., Darbyshire, F., Evans, R.L., Grutter, H., Jones, A.G. & Yuan, X., Mitrovica, J.X. & Forte, A.M., 2004. A new inference of mantle viscosity 2009. The elusive lithosphere–asthenosphere boundary (LAB) beneath based upon joint inversion of convection and glacial isostatic adjustment cratons, Lithos, 109, 1–22. data, Earth planet. Sci. Lett., 225, 177–189. Farrell, W.E. & Clark, J.A., 1976. On postglacial sea level, Geophys. J. R. Mitrovica, J.X. & Milne, G.A., 2003. On post-glacial sea level: I. General astr. Soc., 46, 647–667. theory, Geophys. J. Int., 154, 253–267. Gomez, N., Latychev, K. & Pollard, D., 2018. A coupled ice sheet-sea level Mitrovica, J.X., Wahr, J., Matsuyama, I. & Paulson, A., 2005. The rotational model incorporating 3D Earth structure: variations in Antarctica during stability of an ice-age earth, Geophys. J. Int., 161, 491–506. the last deglacial retreat, J. Clim., doi:10.1175/JCLI-D-17-0352.1. Morelli, A. & Danesi, S., 2004. Seismological imaging of the Antarctic Gomez, N., Mitrovica, J.X., Huybers, P. & Clark, P.U., 2010. Sea level as continental lithosphere: a review, Glob. Planet. Change, 42, 155–165. a stabilizing factor for marine-ice-sheet grounding lines, Nat. Geosci, 3, Nield, G.A. et al., 2014. Rapid bedrock uplift in the Antarctic Peninsula 850–853. explained by viscoelastic response to recent ice unloading, Earth planet. Gomez, N., Pollard, D. & Holland, D., 2015. Sea-level feedback lowers Sci. Lett., 397, 32–41. projections of future Antarctic Ice-Sheet mass loss, Nat. Commun., 6, Nield, G.A., Whitehouse, P.L., King, M.A. & Clarke, P.J., 2016. Glacial 8798, doi:10.1038/ncomms9798. isostatic adjustment in response to changing Late Holocene behaviour of Gomez, N., Pollard, D. & Mitrovica, J.X., 2013. A 3-D coupled ice sheet – ice streams on the Siple Coast, West Antarctica, Geophys. J. Int., 205, sea level model applied to Antarctica through the last 40 ky, Earth planet. 1–21. Sci. Lett., 384, 88–99. Nield, G.A., Whitehouse, P.L., King, M.A., Clarke, P.J. & Bentley, M.J., Grand, S.P., 2002. Mantle shear–wave tomography and the fate of subducted 2012. Increased ice loading in the Antarctic Peninsula since the 1850s slabs, Phil. Trans. R. Soc. Lond., A, 360, 2475–2491. Downloaded from https://academic.oup.com/gji/article/214/2/811/4980921 by DeepDyve user on 12 July 2022 824 G. A. Nield et al. and its effect on glacial isostatic adjustment, Geophys. Res. Lett., 39, van der Wal, W., Whitehouse, P.L. & Schrama, E.J.O., 2015. Effect of GIA L17504. models with 3D composite mantle viscosity on GRACE mass balance Peltier, W.R., 1974. The impulse response of a Maxwell Earth, Rev. Geophys., estimates for Antarctica, Earth planet. Sci. Lett., 414, 134–143. 12, 649–669. Watts, A.B., Zhong, S.J. & Hunter, J., 2013. The behavior of the lithosphere Peltier, W.R., 2004. Global glacial isostasy and the surface of the ice-age on seismic to geologic timescales, Annu. Rev. Earth planet. Sci., 41, 443– Earth: the ICE-5G (VM2) Model and GRACE, Annu. Rev. Earth planet. 468. Sci., 32, 111–149. Wessel, P. & Smith, W.H.F., 1998. New, improved version of generic mapping Peltier, W.R., Argus, D.F. & Drummond, R., 2015. Space geodesy constrains tools released, EOS, Trans. Am. geophys. Un., 79, 579–579. ice age terminal deglaciation: the global ICE-6G C (VM5a) model, J. Whitehouse, P., Latychev, K., Milne, G.A., Mitrovica, J.X. & Kendall, R., geophys. Res., 120, 450–487. 2006. Impact of 3-D Earth structure on Fennoscandian glacial isostatic ad- Pollard, D., Gomez, N. & Deconto, R.M., 2017. Variations of the Antarctic justment: Implications for space-geodetic estimates of present-day crustal ice sheet in a coupled ice sheet-Earth-Sea level model: sensitivity to deformations, Geophys. Res. Lett., 33, doi:10.1029/2006GL026568. viscoelastic Earth properties, J. geophys. Res., 122, 2124–2138. Whitehouse, P.L., Bentley, M.J. & Le Brocq, A.M., 2012a. A deglacial Priestley, K. & McKenzie, D., 2013. The relationship between shear wave model for Antarctica: geological constraints and glaciological modelling velocity, temperature, attenuation and viscosity in the shallow part of the as a basis for a new model of Antarctic glacial isostatic adjustment, Quat. mantle, Earth planet. Sci. Lett., 381, 78–91. Sci. Rev., 32, 1–24. Rignot, E., Mouginot, J., Morlighem, M., Seroussi, H. & Scheuchl, B., 2014. Whitehouse, P.L., Bentley, M.J., Milne, G.A., King, M.A. & Thomas, I.D., Widespread, rapid grounding line retreat of Pine Island, Thwaites, Smith, 2012b. A new glacial isostatic adjustment model for Antarctica: calibrated and Kohler glaciers, West Antarctica, from 1992 to 2011, Geophys. Res. and tested using observations of relative sea-level change and present-day Lett., 41, 3502–3509. uplift rates, Geophys. J. Int., 190, 1464–1482. Ritsema, J., Deuss, A., van Heijst, H.J. & Woodhouse, J.H., 2011. S40RTS: Wolstencroft, M. et al., 2015. Uplift rates from a new high-density GPS a degree-40 shear-velocity model for the mantle from new Rayleigh wave network in Palmer Land indicate signiﬁcant late Holocene ice loss in the dispersion, teleseismic traveltime and normal-mode splitting function southwestern Weddell Sea, Geophys. J. Int., 203, 737–754. measurements, Geophys. J. Int., 184, 1223–1236. Wolstencroft, M., Shen, Z., Tor ¨ nqvist, T.E., Milne, G.A. & Kulp, M., 2014. Sasgen, I. et al., 2017. Joint inversion estimate of regional glacial isostatic Understanding subsidence in the Mississippi Delta region due to sediment, adjustment in Antarctica considering a lateral varying Earth structure ice, and ocean loading: Insights from geophysical modeling, J. geophys. (ESA STSE Project REGINA), Geophys. J. Int., 211, 1534–1553. Res., 119, 3838–3856. Schoof, C., 2007. Ice sheet grounding line dynamics: steady Wu, P., 1999. Modelling postglacial sea levels with power-law rheology and states, stability, and hysteresis, J. geophys. Res., 112, F03S28, a realistic ice model in the absence of ambient tectonic stress, Geophys. doi:10.1029/2006JF000664. J. Int., 139, 691–702. Schotman, H.H.A., Wu, P. & Vermeersen, L.L.A., 2008. Regional pertur- Wu, P., 2004. Using commercial ﬁnite element packages for the study of bations in a global background model of glacial isostasy, Phys. Earth earth deformations, sea levels and the state of stress, Geophys.J.Int., 158, planet. Inter., 171, 323–335. 401–408. Shapiro, N.M. & Ritzwoller, M.H., 2004. Inferring surface heat ﬂux dis- Wu, P. & Johnston, P., 1998. Validity of using ﬂat-earth ﬁnite element models tributions guided by a global seismic model: particular application to in the study of postglacial rebound, in Dynamics of the Ice Age Earth: a Antarctica, Earth planet. Sci. Lett., 223, 213–224. Modern Perspective, pp. 191–202, ed. Wu, P. Trans Tech Pub. Steffen, H., Kaufmann, G. & Lampe, R., 2014. Lithosphere and upper- Zhao, C., King, M.A., Watson, C.S., Barletta, V.R., Bordoni, A., Dell, M. mantle structure of the southern Baltic Sea estimated from modelling & Whitehouse, P.L., 2017. Rapid ice unloading in the Fleming Glacier relative sea-level data with glacial isostatic adjustment, Solid Earth, 5, region, southern Antarctic Peninsula, and its effect on bedrock uplift rates, 447–459. Earth planet. Sci. Lett., 473, 164–176. Steffen, H., Kaufmann, G. & Wu, P., 2006. Three-dimensional ﬁnite-element Zhao, S., Lambeck, K. & Lidberg, M., 2012. Lithosphere thickness and man- modeling of the glacial isostatic adjustment in Fennoscandia, Earth tle viscosity inverted from GPS-derived deformation rates in Fennoscan- planet. Sci. Lett., 250, 358–375. dia, Geophys. J. Int., 190, 278–292. van der Wal, W., Barnhoorn, A., Stocchi, P., Gradmann, S., Wu, P., Drury, Zhong, S., Paulson, A. & Wahr, J., 2003. Three-dimensional ﬁnite-element M. & Vermeersen, B., 2013. Glacial isostatic adjustment model with modelling of Earth’s viscoelastic deformation: effects of lateral variations composite 3-D Earth rheology for Fennoscandia, Geophys. J. Int., 194, in lithospheric thickness, Geophys. J. Int., 155, 679–695. 61–77.

Journal

Geophysical Journal International – Oxford University Press

Published: Aug 1, 2018

References

Computations of the viscoelastic response of a 3-D compressible Earth to surface loading: an application to Glacial Isostatic Adjustment in Antarctica and Canada

A, G; Wahr, J.; Zhong, S.
Temperature, lithosphere–asthenosphere boundary, and heat flux beneath the Antarctic Plate inferred from seismic velocities
S-velocity model and inferred Moho topography beneath the antarctic plate from rayleigh waves
The Antarctica component of postglacial rebound model ICE-6G_C (VM5a) based on GPS positioning, exposure age dating of ice thicknesses, and relative sea level histories

Argus, D.F.; Peltier, W.R.; Drummond, R.; Moore, A.W.
Iceland rising: solid Earth response to ice retreat inferred from satellite radar interferometry and visocelastic modeling

Auriac, A.; Spaans, K.H.; Sigmundsson, F.; Hooper, A.; Schmidt, P.; Lund, B.
Lateral, radial, and temporal variations in upper mantle viscosity and rheology under Scandinavia

Barnhoorn, A.; van der Wal, W.; Vermeersen, B.L.A.; Drury, M.R.
How to evaluate model-derived deglaciation chronologies: a case study using Antarctica

Briggs, R.D.; Tarasov, L.
Variability in the mass flux of the Ross ice streams, West Antarctica, over the last millennium

Catania, G.; Hulbe, C.; Conway, H.; Scambos, T.A.; Raymond, C.F.
Preliminary reference Earth model

Dziewonski, A.M.; Anderson, D.L.
The elusive lithosphere–asthenosphere boundary (LAB) beneath cratons

Eaton, D.W.; Darbyshire, F.; Evans, R.L.; Grütter, H.; Jones, A.G.; Yuan, X.
On

Farrell, W.E.; Clark, J.A.
A coupled ice sheet-sea level model incorporating 3D Earth structure: variations in Antarctica during the last deglacial retreat

Gomez, N.; Latychev, K.; Pollard, D.
Sea level as a stabilizing factor for marine-ice-sheet grounding lines

Gomez, N.; Mitrovica, J.X.; Huybers, P.; Clark, P.U.
Sea-level feedback lowers projections of future Antarctic Ice-Sheet mass loss

Gomez, N.; Pollard, D.; Holland, D.
A 3-D coupled ice sheet – sea level model applied to Antarctica through the last 40 ky

Gomez, N.; Pollard, D.; Mitrovica, J.X.
Mantle shear–wave tomography and the fate of subducted slabs
Upper mantle structure of central and West Antarctica from array analysis of Rayleigh wave phase velocities
Getting Started with ABAQUS, Version (6.14)

Hibbitt, D.; Karlsson, B.; Sorensen, P.
Inside the Subduction Factory

Hirth, G.; Kohlstedt, D.; Eiler, J.
Antarctic contribution to sea level rise observed by GRACE with improved GIA correction

Ivins, E.R.; James, T.S.; Wahr, J.; Schrama, E.J.O.; Landerer, F.W.; Simon, K.M.
On lateral viscosity contrast in the mantle and the rheology of low-frequency geodynamics

Ivins, E.R.; Sammis, C.G.
Lateral viscosity variations beneath Antarctica and their implications on regional rebound motions and seismotectonics

Kaufmann, G.; Wu, P.; Ivins, E.R.
Glacial isostatic adjustment in Fennoscandia for a laterally heterogeneous earth

Kaufmann, G.; Wu, P.; Li, G.
On post-glacial sea level – II. Numerical formulation and comparative results on spherically symmetric models

Kendall, R.A.; Mitrovica, J.X.; Milne, G.A.
Geodetic measurements reveal similarities between post–last glacial maximum and present-day mass loss from the Greenland ice sheet
Lower satellite-gravimetry estimates of Antarctic sea-level contribution

King, M.A.; Bingham, R.J.; Moore, P.; Whitehouse, P.L.; Bentley, M.J.; Milne, G.A.
The influence of viscosity structure in the lithosphere on predictions from models of glacial isostatic adjustment

Kuchar, J.; Milne, G.A.
Influence of lithospheric thickness variations on 3-D crustal velocities due to glacial isostatic adjustment

Latychev, K.; Mitrovica, J.X.; Tamisiea, M.E.; Tromp, J.; Moucha, R.
Glacial isostatic adjustment on 3-D Earth models: a finite-volume formulation

Latychev, K.; Mitrovica, J.X.; Tromp, J.; Tamisiea, M.E.; Komatitsch, D.; Christara, C.C.
Inferences of mantle viscosity based on ice age data sets: Radial structure

Lau, H.C.P.; Mitrovica, J.X.; Austermann, J.; Crawford, O.; Al-Attar, D.; Latychev, K.
An assessment of forward and inverse GIA solutions for Antarctica

Martín-Español, A.; King, M.A.; Zammit-Mangion, A.; Andrews, S.B.; Moore, P.; Bamber, J.L.
Inverting the Fennoscandian relaxation-time spectrum in terms of an axisymmetric viscosity distribution with a lithospheric root

Martinec, Z.; Wolf, D.
Postglacial sea-level change on a rotating Earth: first results from a gravitationally self-consistent sea-level equation

Milne, G.A.; Mitrovica, J.X.
A new inference of mantle viscosity based upon joint inversion of convection and glacial isostatic adjustment data

Mitrovica, J.X.; Forte, A.M.
On post-glacial sea level: I. General theory

Mitrovica, J.X.; Milne, G.A.
The rotational stability of an ice-age earth

Mitrovica, J.X.; Wahr, J.; Matsuyama, I.; Paulson, A.
Seismological imaging of the Antarctic continental lithosphere: a review

Morelli, A.; Danesi, S.
Rapid bedrock uplift in the Antarctic Peninsula explained by viscoelastic response to recent ice unloading
Glacial isostatic adjustment in response to changing Late Holocene behaviour of ice streams on the Siple Coast, West Antarctica

Nield, G.A.; Whitehouse, P.L.; King, M.A.; Clarke, P.J.
Increased ice loading in the Antarctic Peninsula since the 1850s and its effect on glacial isostatic adjustment

Nield, G.A.; Whitehouse, P.L.; King, M.A.; Clarke, P.J.; Bentley, M.J.
The impulse response of a Maxwell Earth
Global glacial isostasy and the surface of the ice-age Earth: the ICE-5G (VM2) Model and GRACE
Space geodesy constrains ice age terminal deglaciation: the global ICE-6G_C (VM5a) model

Peltier, W.R.; Argus, D.F.; Drummond, R.
Variations of the Antarctic ice sheet in a coupled ice sheet-Earth-Sea level model: sensitivity to viscoelastic Earth properties

Pollard, D.; Gomez, N.; Deconto, R.M.
The relationship between shear wave velocity, temperature, attenuation and viscosity in the shallow part of the mantle

Priestley, K.; McKenzie, D.
Widespread, rapid grounding line retreat of Pine Island, Thwaites, Smith, and Kohler glaciers, West Antarctica, from 1992 to 2011

Rignot, E.; Mouginot, J.; Morlighem, M.; Seroussi, H.; Scheuchl, B.
S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements

Ritsema, J.; Deuss, A.; van Heijst, H.J.; Woodhouse, J.H.
Joint inversion estimate of regional glacial isostatic adjustment in Antarctica considering a lateral varying Earth structure (ESA STSE Project REGINA)
Ice sheet grounding line dynamics: steady states, stability, and hysteresis
Regional perturbations in a global background model of glacial isostasy

Schotman, H.H.A.; Wu, P.; Vermeersen, L.L.A.
Inferring surface heat flux distributions guided by a global seismic model: particular application to Antarctica

Shapiro, N.M.; Ritzwoller, M.H.
Lithosphere and upper-mantle structure of the southern Baltic Sea estimated from modelling relative sea-level data with glacial isostatic adjustment

Steffen, H.; Kaufmann, G.; Lampe, R.
Three-dimensional finite-element modeling of the glacial isostatic adjustment in Fennoscandia

Steffen, H.; Kaufmann, G.; Wu, P.
Glacial isostatic adjustment model with composite 3-D Earth rheology for Fennoscandia

van der Wal, W.; Barnhoorn, A.; Stocchi, P.; Gradmann, S.; Wu, P.; Drury, M.; Vermeersen, B.
Effect of GIA models with 3D composite mantle viscosity on GRACE mass balance estimates for Antarctica

van der Wal, W.; Whitehouse, P.L.; Schrama, E.J.O.
The behavior of the lithosphere on seismic to geologic timescales

Watts, A.B.; Zhong, S.J.; Hunter, J.
New, improved version of generic mapping tools released

Wessel, P.; Smith, W.H.F.
Impact of 3-D Earth structure on Fennoscandian glacial isostatic adjustment: Implications for space-geodetic estimates of present-day crustal deformations

Whitehouse, P.; Latychev, K.; Milne, G.A.; Mitrovica, J.X.; Kendall, R.
A deglacial model for Antarctica: geological constraints and glaciological modelling as a basis for a new model of Antarctic glacial isostatic adjustment

Whitehouse, P.L.; Bentley, M.J.; Le Brocq, A.M.
A new glacial isostatic adjustment model for Antarctica: calibrated and tested using observations of relative sea-level change and present-day uplift rates

Whitehouse, P.L.; Bentley, M.J.; Milne, G.A.; King, M.A.; Thomas, I.D.
Uplift rates from a new high-density GPS network in Palmer Land indicate significant late Holocene ice loss in the southwestern Weddell Sea
Understanding subsidence in the Mississippi Delta region due to sediment, ice, and ocean loading: Insights from geophysical modeling

Wolstencroft, M.; Shen, Z.; Törnqvist, T.E.; Milne, G.A.; Kulp, M.
Modelling postglacial sea levels with power-law rheology and a realistic ice model in the absence of ambient tectonic stress
Using commercial finite element packages for the study of earth deformations, sea levels and the state of stress
Dynamics of the Ice Age Earth: a Modern Perspective

Wu, P.; Johnston, P.; Wu, P.
Rapid ice unloading in the Fleming Glacier region, southern Antarctic Peninsula, and its effect on bedrock uplift rates

Zhao, C.; King, M.A.; Watson, C.S.; Barletta, V.R.; Bordoni, A.; Dell, M.; Whitehouse, P.L.
Lithosphere thickness and mantle viscosity inverted from GPS-derived deformation rates in Fennoscandia

Zhao, S.; Lambeck, K.; Lidberg, M.
Three-dimensional finite-element modelling of Earth's viscoelastic deformation: effects of lateral variations in lithospheric thickness

Zhong, S.; Paulson, A.; Wahr, J.

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

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

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

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

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

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

Abstract

Journal

Recommended Articles

References

Our policy towards the use of cookies