## Abstract

The West Antarctic Ice Sheet (WAIS) overlies a thin, variable-thickness lithosphere and a shallow upper-mantle region of laterally varying and, in some regions, very low (~10^{18} Pa s) viscosity. We explore the extent to which viscous effects may affect predictions of present-day geoid and crustal deformation rates resulting from Antarctic ice mass flux over the last quarter century and project these calculations into the next half century, using viscoelastic Earth models of varying complexity. Peak deformation rates at the end of a 25-yr simulation predicted with an elastic model underestimate analogous predictions that are based on a 3D viscoelastic Earth model (with minimum viscosity below West Antarctica of 10^{18} Pa s) by ~15 and ~3 mm yr^{−1} in the vertical and horizontal directions, respectively, at sites overlying low-viscosity mantle and close to high rates of ice mass flux. The discrepancy in uplift rate can be reduced by adopting 1D Earth models tuned to the regional average viscosity profile beneath West Antarctica. In the case of horizontal crustal rates, adopting 1D regional viscosity models is no more accurate in recovering predictions that are based on 3D viscosity models than calculations that assume a purely elastic Earth. The magnitude and relative contribution of viscous relaxation to crustal deformation rates will likely increase significantly in the next several decades, and the adoption of 3D viscoelastic Earth models in analyses of geodetic datasets [e.g., Global Navigation Satellite System (GNSS); Gravity Recovery and Climate Experiment (GRACE)] will be required to accurately estimate the magnitude of Antarctic modern ice mass flux in the progressively warming world.

## 1. Introduction

Climate is changing, with warming that leads to an increase in ice melting and corresponding global mean sea level rise. Constraining and projecting sea level variability requires an accurate estimate of the size and geometry of the meltwater sources. One way to quantify ice mass flux is to measure the associated deformation of Earth. When ice melting occurs, the resulting (ice plus ocean) mass redistribution perturbs Earth’s gravitational field and solid surface, and these effects can be measured using a suite of geodetic methods, including, for example, satellite gravity observations and surveying using the Global Navigation Satellite System (GNSS). Such analyses require assumptions of Earth’s rheological response time to the loading. On shorter time scales (decades–centuries) Earth’s response is often assumed to be primarily elastic, whereas on longer (e.g., ice age) time scales, Earth’s response is clearly viscoelastic. In this study, we use the West Antarctic as a case study of the transition between these two regimes. In particular, we adopt a series of ice-melting scenarios extending over the past 25 years and projecting into the next half century and incorporate 3D viscoelastic Earth structure in glacial isostatic adjustment (GIA) simulations to explore the time scale over which viscous forces become significant in driving gravity perturbations and 3D crustal motions in the South Pole region.

Geodetic measurements can provide important constraints on ice melt. For example, the Gravity Recovery and Climate Experiment (GRACE), including a satellite mission operational from March 2002 to October 2017 and the current follow-on mission GRACE-FO (launched May 2018), maps geoid anomalies into surface mass changes assuming an elastic response of the solid Earth (Wahr et al. 1998). An analysis of GRACE data extending over the period April 2002 to January 2009 has argued that the West Antarctic lost an average of 132 ± 26 Gt of ice per year during this period (Chen et al. 2009). This ice melt signal, particularly in the Amundsen Sea sector, appears to have accelerated since 2002 (e.g., Velicogna et al. 2014; Shepherd et al. 2018). However, substantial uncertainty in these estimates comes from the contribution to mass changes from movement of the solid Earth associated with past loading over the last glacial cycle (i.e., GIA). GNSS measurements in Antarctica are also sensitive to modern (annual–century time scale) ice mass change and various approaches have been used to separate this signal from the GIA component of the crustal deformation. Within these analyses, the crustal response to modern melting has been computed using either purely elastic Earth models (e.g., Bevis et al. 2009; Thomas et al. 2011; Argus et al. 2014; Martín-Español et al. 2016; Caron et al. 2018; Schumacher et al. 2018) or by augmenting these calculations to include viscous relaxation to explain anomalously rapid uplift rates in specific areas of the West Antarctic (e.g., Nield et al. 2014; Zhao et al. 2017; Barletta et al. 2018). With few exceptions (e.g., Argus et al. 2014; Zhao et al. 2017), only the vertical component of GNSS measurements has been considered in such analyses.

Antarctica has a complex geologic setting that is not well described by purely elastic models. The East Antarctic is characterized by an old, cold craton with a thick lithosphere in excess of 200 km (Morelli and Danesi 2004; Heeszel et al. 2013), while the West Antarctic is dominated by a failed rift system (Wörner 1999) that has thinned the lithosphere to <100 km (An et al. 2015a; Heeszel et al. 2016). Simple, thermal interpretations of seismic tomographic images of the region (Ritzwoller et al. 2001; Morelli and Danesi 2004; Hansen et al. 2014; Lloyd et al. 2015; Heeszel et al. 2016) suggest that mantle viscosities below parts of West Antarctica are also significantly lower than both the regional Antarctic and global average, with viscosities suggested to be as low as ~10^{18} Pa s under, for example, volcanic Marie Byrd Land (Kaufmann et al. 2005, Hay et al. 2017). While uncertainty remains in mapping seismic wave speed anomalies to viscosity structure, these estimates are consistent with inferences of low asthenospheric-mantle viscosity based on analyses of GNSS-determined crustal uplift rates in the Antarctic Peninsula (Nield et al. 2014; Zhao et al. 2017) and Amundsen Sea Embayment region (Barletta et al. 2018). Analysis of xenoliths collected from Marie Byrd Land suggest that local viscosities in the shallow mantle below this area may be as low as 10^{16} Pa s (Chatzaras et al. 2016; S. C. Kruckenberg 2019, personal communication).

This complicated 3D structure has already been studied in the context of GIA in response to the last ice age (Kaufmann et al. 2005; A et al. 2013; van der Wal et al. 2015; Gomez et al. 2018). However, such low viscosities indicate Maxwell times of less than a year, suggesting that viscous effects play a role even in the response to modern melting over the West Antarctic. Previous studies have considered the impact of viscous relaxation on the response to modern melting at other sites characterized by shallow mantle viscosities of order 10^{18} Pa s. These include examinations of the crustal response to melting in the Antarctic Peninsula (Nield et al. 2014; Zhao et al. 2017), the Amundsen Sea Embayment (Barletta et al. 2018), Patagonia (Richter et al. 2016), Iceland (Auriac et al. 2013), eastern Greenland (Khan et al. 2016), and Alaska (Tamisiea et al. 2005; James et al. 2009). This issue also motivated the study of Hay et al. (2017), who explored the impact of 3D low-viscosity structure beneath the West Antarctic on the sea level fingerprints of modern melting and concluded that peak sea level fall in the West Antarctic associated with a local melting event of duration 25 years will increase by 25% relative to a purely elastic simulation.

We have two goals in this study. First, we seek to estimate the contribution of 3D variations in mantle viscosity beneath West Antarctica to predictions of the gravitational field and crustal deformation response to ice mass flux over the past 25 years and projected forward over the next half century. Second, given the significant technical requirements involved in treating 3D viscoelastic Earth structure in such loading calculations, we explore whether 1D viscosity models can be found that provide a reasonable approximation of these 3D effects.

## 2. Methods

To predict Earth’s response to a change in ice loading, one must account, in a gravitationally self-consistent manner, for the flux of water into and out of the ocean basins. In the present study we adopt the sea level theory described by Gomez et al. (2010). This theory assumes the initial topography is known, and it incorporates effects associated with time-varying shoreline geometry (Johnston 1993; Milne and Mitrovica 1998; Mitrovica and Milne 2003) and load-induced Earth rotation variations (Mitrovica et al. 2005).

In studies adopting 1D viscoelastic Earth models, loading calculations inherent to the sea level theory are usually based on viscoelastic Love number theory (Peltier 1974). The incorporation of 3D Earth structure requires a more complex treatment of load-induced perturbations to the gravitational field and crust, and in this regard we adopt the finite-volume treatment of Latychev et al. (2005). With recent improvements (e.g., Hay et al. 2017; Gomez et al. 2018), we extend the treatment to include a laterally varying resolution in the computational grid to accommodate available regional models of higher spatial resolution. The global model we adopt is characterized by an average spatial (horizontal and vertical) resolution of 12 km to the base of the crust, 25 km to a depth of 220 km, and 50 km to the core mantle boundary. The regional model, which asymmetrically covers the Antarctic plate spatially, extends to depths of approximately 350 km, and it is characterized by an average spatial (horizontal and vertical) resolution of 5 km to the base of the crust, 12 km to a depth of 220 km, and 25 km to a depth of 350 km (see Fig. S1 in the online supplemental material). We use the finite-volume software and gridding scheme in all calculations presented in this study, including those in which Earth structure varies only with depth.

Our calculations require two inputs: a model for Earth’s viscoelastic structure, and the space–time history of ice cover. We describe these inputs below.

### a. Earth models

We consider a suite of Earth models in this study. All Earth models assume a Maxwell viscoelastic mantle rheology that is compressible in the elastic limit. The elastic and density structures of the models are provided by the 1D seismic Preliminary Reference Earth Model (PREM; Dziewonski and Anderson 1981).

The first Earth model is purely elastic (*M*_{EL}). In this case, the computed perturbations to the solid surface and gravitational field do not depend on the duration of the simulation, only on the net change in ice volume between the beginning and end of the calculation. The second model (*M*_{1D}) has viscoelastic structure that varies with depth alone. In particular, the 96-km-thick lithosphere overlies uniform upper and lower mantle viscosities of 5 × 10^{20} and 5 × 10^{21} Pa s, respectively. This viscosity profile is shown in Fig. 1a (black curve). The *M*_{1D} model is characteristic of 1D viscosity models favored in most GIA-based inferences of mantle viscosity based on globally distributed datasets (e.g., Lambeck et al. 1998; Mitrovica and Forte 2004; Lau et al. 2016).

The third model (*M*_{3D}) is defined by an elastic lithosphere of variable thickness (Fig. 1b) and 3D mantle viscosity structure (Fig. 1c). Globally, we adopt the spatially varying lithospheric thickness model of Conrad and Lithgow-Bertelloni (2006), but within the Antarctic plate we use the higher-resolution lithospheric model of An et al. (2015a). The full model is scaled to yield a global mean lithospheric thickness of 96 km. The 3D mantle viscosity structure of *M*_{3D} is built from three different seismic tomography studies that span global to regional (Antarctic) scale (Ritsema et al. 2011; Heeszel et al. 2016; An et al. 2015b). The model, which is described in full detail in Hay et al. (2017), involves a free parameter that controls the level of lateral variability in mantle viscosity. In our standard run, this parameter is chosen such that the Earth model is characterized by a five-order-of-magnitude (peak to peak) variation in viscosity in the asthenosphere beneath East and West Antarctica, where the latter region has a minimum viscosity of ~10^{18} Pa s. To test the sensitivity of the results to this choice, we consider two other values of the free parameter that yield minimum viscosities of 10^{17} and 10^{19} Pa s beneath West Antarctica (models *M*_{3D-L} and *M*_{3D-H}, respectively). Globally, the 3D Earth models are constrained to have a spherically averaged depth profile that matches the model *M*_{1D}. The results are relatively insensitive to this choice of spherically averaged viscosity given that we tune our model free parameter to yield specific lower bounds on viscosity (10^{17}, 10^{18}, or 10^{19} Pa s) beneath the West Antarctic.

The fourth model (*M*_{MBL}) is a second 1D viscoelastic model constructed from the regional mantle viscosity structure beneath Marie Byrd Land in model *M*_{3D}, the area overlying the asthenospheric low-viscosity zone in the *M*_{3D} model. The model has a 71-km-thick lithosphere and a highly variable mantle viscosity profile (Fig. 1a, red curve) constructed by taking the cylindrical average of viscosity with depth from the surface to the core mantle boundary, using a 555-km-radius circle centered on 79°S, 124°W (Fig. 1c).

### b. Ice model

The normalized, uniform melting, or full collapse scenario generally used to calculate sea level fingerprints, including in the study of Hay et al. (2017), does not capture the complex geometry or magnitude of recent mass loss in Antarctica. Ice mass redistribution occurs primarily via ice streams within the ice sheet and via calving and melting at the periphery of the ice sheet (Bennett 2003; Shepherd et al. 2018). Antarctica contributed 0.27 ± 0.11 mm yr^{−1} equivalent sea level rise over the period from 1993 to 2010 (Vaughan et al. 2013) and this rate has increased considerably since 2010 (Harig and Simons 2015; Martín-Español et al. 2016; Shepherd et al. 2018).

In the present study, we adopt a suite of ice melt histories over the Antarctic. The first two models extend over 25 years (1992–2017), consistent with the period over which the modern Antarctic Ice Sheet has been significantly out of mass balance, that is, from 1992 to present day (Shepherd et al. 2018). The first ice history (I-GR) is based on geographically variable melt rates inferred from GRACE satellite gravity data collected from 2003 to 2014 (Harig and Simons 2015). Specifically, the geometry of the ice melt is defined by the mean annual change in ice thickness over that time period (Fig. 2a), and we apply a constant melt volume change equivalent to 0.26 mm yr^{−1} of global mean sea level rise over the entire 25-yr simulation (Fig. 2c, inset; Harig and Simons 2015).

Our second ice history (I-ME) adopts the full spatiotemporal evolution of the Martín-Español et al. (2016) reconstruction from 2003 to 2013. (Fig. 2b shows the mean annual ice thickness change across the entire period.) For the period 1992–2002 we use the 2003 mass flux geometry in the Martín-Español et al. (2016) reconstruction and we follow the integrated, time-varying mass flux inferred by the IMBIE team for this 10-yr period (Shepherd et al. 2018; Fig. 2c, inset). For 2014–17, the I-ME ice history adopts the 2013 mass flux geometry in the Martín-Español et al. (2016) reconstruction and follows the integrated mass flux inferred by Shepherd et al. (2018) for the same 4-yr period (Fig. 2c, inset).

Our third ice history (I-MEG) is identical to I-ME over the period 1992–2017, but it extends this history for 50 additional years, to 2067. Over this latter period, the geometry of the mass flux is held fixed, and the integrated magnitude of the flux follows the trend predicted by Golledge et al. (2019) in their coupled ice sheet/ice shelf simulations of the Antarctic Ice Sheet over the twenty-first century (Fig. 2c).

## 3. Results and discussion

The Maxwell time associated with the model *M*_{1D} is of the order of a millennium and, as a consequence, all of the 25-yr simulations we performed yielded negligible differences between predictions based on the *M*_{1D} and *M*_{EL} Earth models. We therefore omit results generated using the *M*_{1D} model in the figures discussed below. For the purposes of this study, the *M*_{1D} model, a standard mantle viscosity profile inferred from GIA datasets, is essentially indistinguishable from a purely elastic Earth model.

In the following sections, we plot predictions based on the model *M*_{3D} and the difference in predictions based on the pair of models (*M*_{3D}, *M*_{EL}) and (*M*_{3D}, *M*_{MBL}). The first of these pairs represents the viscous signal embedded within the 3D Earth model simulation. The difference in the second pair quantifies the extent to which the regional, 1D viscosity model *M*_{MBL} captures the viscous signal within the *M*_{3D}-based simulation.

### a. Geoid rate predictions

Figure 3a shows the rate of change in the height of the geoid at the end of the 25-yr simulation computed using the *M*_{3D} model and the I-GR ice history. This rate, which incorporates perturbations associated with both the ice mass flux and the associated adjustment of the solid Earth, has a peak negative value of ~3.5 mm yr^{−1} over the West Antarctic and a peak positive value of 0.6 mm yr^{−1} over the East Antarctic. Figures 3b and 3c show the difference in predictions of the geoid height rate change generated using the 3D Earth model prediction and the two 1D models *M*_{EL} and *M*_{MBL}, respectively. The magnitude of the peak geoid rate over the West Antarctic in Fig. 3a is 6% smaller (~0.2 mm yr^{−1}) than the analogous predictions based on the model *M*_{EL}, reflecting an increased compensation of the geoid signal (due to ice mass loss) associated with uplift of the crust due to viscous mantle flow relative to the compensation computed using a purely elastic Earth model. The difference in the peak rate over the West Antarctic predicted using the *M*_{MBL} and *M*_{3D} models is less, ~1% (Fig. 3c), indicating that the 1D model tuned to the regional viscosity variation beneath the West Antarctic accurately captures the 3D Earth model prediction.

The results in Figs. 3d–f are analogous to the top row of the figure, but based on the more spatially resolved ice history I-ME. The peaks in the predicted signal based on the model *M*_{3D} are more localized, reflecting the geometry of the underlying mass flux (e.g., Fig. 2b), and the amplitudes are significantly higher. The viscous signal in the geoid rate peaks at 0.37 mm yr^{−1} (Fig. 3e; 1.3% of the peak in Fig. 3d, and, as in the predictions based on the ice history I-GR, this signal is captured to within 0.3% accuracy with the 1D viscoelastic model *M*_{MBL} derived from regional structure below the West Antarctic (Fig. 3f).

These results indicate that low-viscosity structure beneath the West Antarctic has a relatively small impact on predictions of geoid rate, and that analyses of GIA-corrected GRACE measurements over the West Antarctic do not incur significant errors in assuming that modern ice mass loss drives a purely elastic solid Earth response.

### b. Crustal deformation rate predictions

Next, we consider predictions of crustal deformation rates computed using ice history I-ME (top rows of Figs. 4 and 5, and all of Fig. 6). The vertical component of these rates has served as the primary dataset in analyses of GNSS measurements across the Antarctic.

Crustal uplift rates predicted using the 3D Earth model *M*_{3D} (Fig. 4a) are characterized by a peak value of 90 mm yr^{−1} at the location of greatest ice mass flux in the I-ME history, near the Amundsen Sea Embayment (Fig. 2b). The viscous component of the uplift field peaks at 14 mm yr^{−1} (Fig. 4b). The ratio of Figs. 4b and 4a indicates that the viscous signal reaches ~20% of the full calculation in regions where significant uplift rates are predicted (Fig. 5a). A significant component of this viscous signal is captured in the calculation based on the 1D, regionally inferred viscoelastic Earth model *M*_{MBL}; within the zone of pronounced ice melting in the West Antarctic, the discrepancy between predictions based on models *M*_{MBL} and *M*_{3D} (Fig. 4c) ranges from −5.6 to 3.3 mm yr^{−1} (cf. Figs. 5a and 5b; we return to this point below).

These results demonstrate that adopting elastic Earth models to correct GNSS measurements of vertical crustal rates for the signal due to modern melting in the West Antarctic will underestimate the magnitude of the correction, and thus overestimate the residual contribution from the GIA process. Alternatively, using a GIA-corrected field of measured crustal uplift to estimate modern mass flux would overestimate this flux if a purely elastic response is adopted to compute the response to recent melting.

Next, we turn to predictions of horizontal crustal motions predicted at the end of the 25-yr I-ME simulations. Figure 4d shows the results based on model *M*_{3D}; Fig. 4e shows the contributions to this field from viscous effects (i.e., the vector difference of predictions based on models *M*_{EL} and *M*_{3D}), while Fig. 4f is the difference in predictions based on models *M*_{MBL} and *M*_{3D}.

Within the zone of crustal uplift (Fig. 4a), the horizontal rate predictions based on the 3D Earth model *M*_{3D} emanate outward from the zone of peak uplift (Fig. 4d), with peak rates that reach 15 mm yr^{−1}. A comparison of this prediction with the viscous signal (Fig. 4e; note the different scale of the arrows in the two panels) indicates that the outward pattern in this region is dominated by elastic flexure (James and Morgan 1990). Nevertheless, the viscous signal, which drives horizontal deformation inward toward the areas of melt (James and Morgan 1990), exceeds 3 mm yr^{−1} within the zone of crustal uplift, and remains above 1 mm yr^{−1} well outside this region, particularly in oceanic crust to the north (Fig. 4e). The prediction based on the 1D regional viscosity model *M*_{MBL} within the West Antarctic fails to capture the viscous signal in the 3D simulation (i.e., the residuals in Fig. 4f are of similar magnitude to the viscous signal in Fig. 4e) and is, in general, comparable in performance to the elastic model as an approximation to the 3D viscoelastic simulation based on *M*_{3D}. This conclusion is further reinforced when considering results within Marie Byrd Land alone (see Fig. S2 and caption in the online supplemental material).

Figure 6 tracks predicted crustal rates at three representative GNSS sites (Figs. 6b–d; locations are shown in Fig. 6a) in the West Antarctic. (Rates are computed using a sliding window of 5 years.) These sites lie on an east–west arc that spans the zone used to average viscosity in the construction of the 1D regional model *M*_{MBL} (cf. the location of the three sites in Fig. 6a and the dashed circle in Fig. 1c). We show the prediction generated using the 3D viscoelastic Earth model *M*_{3D} (top row of each panel; black lines) and the difference in the predictions (bottom row on each panel) based on model *M*_{3D} and either model *M*_{EL} (blue lines) or model *M*_{MBL} (red lines). The figure provides a measure of the progression in the difference over the full 25-yr time window between a prediction based on the 3D viscoelastic model *M*_{3D} and a prediction that either adopts a purely elastic response of the solid Earth or models the viscous response using a 1D regionally inferred viscosity profile.

Site 5 (GNSS site TOMO) lies closest to the region of highest ice mass flux and both the prediction of crustal uplift and the viscous component of this signal are the largest of the three sites; the viscous component reaches 8.8, 0.3, and 0.5 mm yr^{−1} for crustal rates of uplift, eastward, and northward horizontal deformation, respectively, at the end of the 25-yr simulation. The prediction based on the 1D model tuned to the regional viscosity profile beneath this region of the West Antarctic, *M*_{MBL,} performs well in recovering the uplift rate signal generated with the model *M*_{3D}, but the discrepancy between predictions of horizontal rates between these two models (*M*_{MBL} and *M*_{3D}) is significantly larger than the viscous signal (i.e., *M*_{3D} − *M*_{EL}) after 25 years (1.1 vs 0.3 mm yr^{−1} in the eastward direction, and 1.6 vs 0.5 mm yr^{−1} in the northward direction). That is, one would incur a greater error using the regional 1D viscosity model *M*_{MBL} than a purely elastic model in predicting horizontal crustal rates at this site computed using the 3D viscoelastic model *M*_{3D}.

The viscous signal at site 3 (CLRK) is of order ~1 mm yr^{−1} or less, and simulations based on model *M*_{MBL} have more success in capturing this component of the *M*_{3D}-based horizontal crustal response. For example, at the end of 25 years, the viscous signal in the *M*_{3D} response for the three crustal deformation components is 1.1, 0.9, and 0.3 mm yr^{−1}, respectively, while the analogous predictions for the *M*_{MBL} simulation (i.e., *M*_{3D} − *M*_{MBL}) are lower: −0.1, 0.5, and 0.2 mm yr^{−1}.

Site 4 (SDLY) is closest to the mantle region of lowest viscosity in the Earth model *M*_{3D} (see Fig. 1c). As for site 3, the model *M*_{MBL} is able to capture nearly all of the viscous component of crustal uplift at site 4, and a substantial fraction of the component for the horizontal-north rate, associated with the prediction based on the 3D Earth model *M*_{3D}; however, it does only marginally better than the model *M*_{EL} in predicting the horizontal-east rate computed using the 3D viscoelastic model *M*_{3D}.

Clearly, the magnitude of the viscous response in predictions of 3D crustal rates, and the ability of the 1D model *M*_{MBL} to recover this response, will depend on the location of the site relative to both the geometry of the ice mass flux and the detailed variability in viscosity and lithospheric thickness that characterizes any 3D model prediction. Table 1 explores this issue further by showing predicted 3D crustal rates at all 10 GNSS sites in Fig. 6a at the end of the 25-yr simulation for model *M*_{3D} and the differences *M*_{3D} − *M*_{EL} and *M*_{3D} − *M*_{MBL}.

### c. Sensitivity analysis: Varying the minimum viscosity

Next, we repeat the calculations based on the I-ME ice history, but vary the mapping from seismic velocity anomalies to viscosity to consider Earth models in which the minimum viscosity below the West Antarctic is reduced to 10^{17} Pa s and increased to 10^{19} Pa s (models *M*_{3D-L} and *M*_{3D-H}, respectively). The viscous signal in predictions of the peak magnitude of geoid rate, crustal uplift rate, and the two components of the horizontal crustal rate (i.e., the difference in the peak magnitude of these quantities computed using the set of 3D Earth models and the *M*_{EL} model) as a function of the minimum viscosity below West Antarctica is summarized in Fig. 7. The viscous signal in the peak geoid rate is less than 1 mm yr^{−1} for all three cases. In contrast, the viscous signal in the crustal deformation rates ranges from 6 to 32 mm yr^{−1} for uplift and from 1 to 9 mm yr^{−1} for horizontal deformation.

### d. Sensitivity analysis: Projections across the next 50 years

As a final analysis, we perform a simulation that extends the calculation based on ice history I-ME for an additional 50 years using the global mean sea level (GMSL) trend predicted by Golledge et al. (2019) for the Antarctic Ice Sheet. The Golledge et al. (2019) projection of Antarctic ice mass flux was generated using a coupled ice sheet–ice shelf model forced with a climatology based on CMIP5 outputs, with additional ice sheet–climate feedbacks, and in the period 2017–67 it projects a GMSL rise of ~60 mm (Fig. 2c, main plot). To be consistent with our construction of the I-ME model, we assume that the ice melt geometry across this 50-yr period is given by the mass flux in the final year of the Martín-Español et al. (2016) reconstruction, and we scale the total melt to follow the GMSL curve of Golledge et al. (2019). We denote the model as I-MEG, and emphasize that the first 25 years of the 75-yr ice history are identical to model I-ME.

The bottom two rows of Fig. 5 show results analogous to the top row—predictions of crustal uplift rates at the 25-yr mark of the I-ME simulation—at years 50 and 75 of the I-MEG simulations. Once again, the left panel on each row represents the contribution, in percent, of the viscous signal relative to the signal predicted using the 3D viscoelastic model *M*_{3D} (i.e., *M*_{3D} prediction minus *M*_{EL} prediction, divided by the former). The right panel provides a measure of the ability of the 1D, regionally tuned model, *M*_{MBL}, to capture these viscous effects (i.e., *M*_{3D} prediction minus *M*_{MBL} prediction, divided by the former). In the case of the right column, one should focus on the region close to Marie Byrd Land since the 1D viscosity profile was based on averaging the viscosity below this region (Fig. 1c, dashed circle). However, the large discrepancies evident at other sites in the West Antarctic (Figs. 5d,f) emphasizes that a 1D viscosity model derived from mantle viscosity variations below one region cannot be interpreted as an appropriate model for computing the response in the West Antarctic as a whole.

Moving down the left column of Fig. 5 indicates that viscous effects peak at 20%, 35%, and 55% of the signal in the prediction of crustal uplift rates based on the 3D viscoelastic model *M*_{3D} at years 25 (i.e., calendar year 2017, as discussed above), 50, and 75 of the simulation. A comparison of these values with the results in the right column indicates that using the 1D viscosity model *M*_{MBL} in place of *M*_{EL} captures only about half of this viscous signal near the zone of major ice mass flux.

Figure 8 and Table 2 are analogous to Fig. 6 and Table 1, except that the figure tracks predictions of crustal rates at the same three GNSS sites for the entire 75-yr duration of the ice history I-MEG. The conclusions derived from Fig. 6 for sites within the region of significant mass flux—namely, that the regional model *M*_{MBL} does a reasonable job at capturing the viscous effects in crustal uplift rates predicted using the 3D model *M*_{3D}, and that the same is not in general true for horizontal rates—continue to hold across the longer simulation. We note also that the viscous signals (blue lines in Figs. 8b–d) and the residuals between predictions based on models *M*_{3D} and *M*_{MBL} increase monotonically over time for all three crustal rate components.

## 4. Conclusions

The Antarctic Ice Sheet is a central focus of studies investigating the impact of global climate change, and geodetic measurements, including GRACE satellite gravity data and surveying using GNSS, play a key role in many such studies. These studies follow two distinct approaches. First, the geodetic measurements are corrected for the ongoing influence of the ice age (i.e., GIA) and the residual signal is used to estimate modern ice mass flux. Second, an independent estimate of modern ice mass flux is used to correct the observational data, leaving a signal that is analyzed to constrain the GIA process. In both these approaches, a mapping is required between modern ice mass flux and perturbations to the Earth system associated with this flux.

The goal of the present study has been to use a 3D model of mantle viscosity to quantify the impact of viscous relaxation of the solid Earth within the Antarctic region on predictions of geoid height changes and crustal deformation rates driven by modern melting, a component of the response that has sometimes been neglected in previous work. Our analysis has involved simulations of duration 25 and 75 years; the former is consistent with the period during the modern over which mass flux from the Antarctic is thought to have been significant (Shepherd et al. 2018), and the latter allows us to estimate the viscous signal associated with Antarctica’s projected melting (Schlegel et al. 2018; Golledge et al. 2019; Bulthuis et al. 2019) as Earth moves further into a warming world. Moreover, we have considered a series of ice histories, and quantified the extent to which 1D models of mantle viscosity can accurately account for viscous effects.

We have found that the viscous signal in predictions of peak geoid height changes in a laterally varying Earth model (i.e., *M*_{3D}) are at the level of 0.5 mm yr^{−1} at the end time of the 25-yr simulations, and conclude that studies analyzing existing GRACE gravity data by assuming that modern mass flux drives a purely elastic response of the solid Earth will marginally overestimate the associated geoid signal at this level. This minor level of inaccuracy can be decreased by modeling the geoid height changes using a 1D Earth model with a radial viscosity profile tuned to the regional viscosity variation under the West Antarctic. The ability of a regionally tuned, 1D viscosity profile to accurately reproduce the present-day geoid signal computed using an Earth model with lateral variations in Earth structure has also been demonstrated in studies of glacial isostatic adjustment in response to the last ice age (Paulson et al. 2005).

Our results demonstrate that viscous effects on crustal deformation rates due to modern ice mass flux are significantly larger than the corresponding effects on the geoid height. A et al. (2013) reached the same conclusion for signals associated with glacial isostatic adjustment. The viscous signal in peak uplift rates predicted using 3D viscoelastic Earth models at sites near areas of significant ice melting ranges from 6 to 32 mm yr^{−1} at the 25-yr mark in the simulations when the minimum viscosity beneath the region is varied from 10^{19} to 10^{17} Pa s (Fig. 7b). A portion of this viscous signal can be modeled by adopting a 1D viscosity profile more characteristic of the low-viscosity region beneath the West Antarctic (*M*_{MBL}). Specifically, results from our suite of simulations extending 25 to 75 years indicate that the *M*_{MBL} model captures approximately 50% of the viscous signal predicted using the 3D viscoelastic Earth model *M*_{3D} (Fig. 5).

The inadequacy of 1D Earth models in reproducing horizontal crustal rates computed using a more complex Earth model with 3D variations in viscoelastic structure is pronounced (Figs. 4, 6, and 8), and the discrepancy between the two will depend on the details of the melt geometry and rheological variability. Indeed, at some GNSS sites within zones of high ice mass flux in the West Antarctic, we have found that the 1D viscosity profile derived from a regional averaging of the viscosity field below the region (*M*_{MBL}) yields horizontal crustal rate predictions that are more discrepant from predictions based on the 3D viscoelastic Earth model *M*_{3D} than predictions generated using an elastic model (*M*_{EL}) (Fig. 6d).

We conclude that horizontal crustal motions due to modern ice mass flux in the Antarctic cannot in general be accurately modeled using any 1D Earth model, and thus the signal due to modern melting in GNSS derived horizontal motions in the West Antarctic should be analyzed using models that incorporate the full complexity of viscoelastic Earth structure beneath the region. This conclusion may have relevance to the results of Zhao et al. (2017), who were unable to find a 1D viscoelastic model that satisfactorily fit their horizontal crustal rate observations. The sensitivity of these observations to lateral viscosity structure suggests that they have limited utility in constraining 1D (i.e., depth varying) viscoelastic structure.

At the end of the 25-yr simulation based on the ice history I-ME and the 3D viscoelastic Earth model *M*_{3D}, the viscous signal within the zone of major ice mass flux reached 14 mm yr^{−1} for uplift rates and 2.5 mm yr^{−1} for horizontal rates. These values exceed average uncertainties in GNSS estimates of these motions (e.g., Barletta et al. 2018). This signal is systematic, not random (Figs. 4–8), and this suggests that viscous effects should be incorporated in predictions of the response to the modern melt signal in regions where shallow mantle viscosities are less than ~10^{19} Pa s now that the time scale of melting from the region has reached one-quarter century (or more).

As a final point, the ice history I-ME and the Earth model *M*_{3D} were adopted as representative of a growing suite of spatially high-resolution reconstructions of ice mass flux in the West Antarctic and 3D mantle viscosity field below this region. Our goal was not to present accurate predictions of crustal deformation rates at specific sites within the West Antarctic, but rather to explore the range of viscous contributions to these predictions in a region of such complexity, and also to assess the ability of regionally tuned, 1D viscosity profiles to capture the viscous signal. Any prediction of crustal rates or geoid anomalies in the West Antarctic will be sensitive to the choice of ice history (including the ice age component of this history) and the details of the Earth model, and efforts to improve constraints on either of these inputs based on geodetic observations must address their coupled nonuniqueness. As recent articles have demonstrated (e.g., Nield et al. 2014; Zhao et al. 2017; Barletta et al. 2018), and the present study emphasizes, in the West Antarctic this effort must contend with viscous effects in the response of the solid Earth to melting over the last few decades. [The viscous signal would be even more pronounced if one adopted an Antarctic loading history that extended over the twentieth century, as in Barletta et al. (2018).] This complication will become ever more relevant and important in the future as we continue to use space-based geodetic techniques to monitor the stability of the Antarctic Ice Sheet in a progressively warming world.

## Acknowledgments

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grants DGE1144152 and DGE1745303, NASA under Grant NNX17AE17G, the National Science and Engineering Research Council, the Canada Research Chairs Program, the Canadian Foundation for Innovation, McGill University, Boston College, and Harvard University.

## REFERENCES

*Climate Change 2013: The Physical Science Basis*, T. F. Stocker et al., Eds., Cambridge University Press, 317–382.

## Footnotes

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-19-0479.s1.

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).