Sea level fingerprints associated with rapid melting of the West Antarctic Ice Sheet (WAIS) have generally been computed under the assumption of a purely elastic response of the solid Earth. The authors investigate the impact of viscous effects on these fingerprints by computing gravitationally self-consistent sea level changes that adopt a 3D viscoelastic Earth model in the Antarctic region consistent with available geological and geophysical constraints. In West Antarctica, the model is characterized by a thin (~65 km) elastic lithosphere and sublithospheric viscosities that span three orders of magnitude, reaching values as low as approximately 4 × 1018 Pa s beneath WAIS. Calculations indicate that sea level predictions in the near field of WAIS will depart significantly from elastic fingerprints in as little as a few decades. For example, when viscous effects are included, the peak sea level fall predicted in the vicinity of WAIS during a melt event will increase by about 20% and about 50%, relative to the elastic case, for events of duration 25 and 100 yr, respectively. The results have implications for studies of sea level change due to both ongoing mass loss from WAIS over the next century and future, large-scale collapse of WAIS on centennial-to-millennial time scales.
The geographic pattern of sea level change following the melting of ice sheets and glaciers is a strong function of the location and time scale of the ice mass change (Farrell and Clark 1976). In the special case of relatively rapid mass loss (i.e., melt events spanning decadal to centennial time scales), these geometries have come to be known as sea level fingerprints. While this specific class of variability has been recognized since the late 1800s (Woodward 1888), the so-called fingerprinting of sea level change received renewed attention beginning at the turn of the current century, following efforts to use widely distributed sea level observations to infer meltwater sources in both the modern (e.g., Mitrovica et al. 2001; Plag and Jüttner 2001; Plag 2006) and ice age (Clark et al. 2002) world.
Modern predictions of sea level fingerprints date to the work of Clark and Lingle (1977), who applied the gravitationally self-consistent sea level theory of Farrell and Clark (1976) to predict sea level changes in response to future melting of the West Antarctic Ice Sheet (WAIS). These predictions assumed a nonrotating Earth and fixed shorelines; the same governing sea level theory was adopted in many subsequent studies (e.g., Clark and Primus 1987; Conrad and Hager 1997). The assumption of fixed shorelines neglects migration due to either onlap or offlap of water at shorelines experiencing local sea level variations or changes in the perimeter of grounded, marine-based ice sheets.
Mitrovica et al. (2001) and Tamisiea et al. (2001) included the feedback of rotational effects into their predictions of sea level fingerprints; that is, their sea level theory included a term for the response to a centrifugal driving potential associated with load-induced polar motion. Moreover, Mitrovica et al. (2009) and Gomez et al. (2010) adopted a generalized sea level theory to incorporate shoreline migration into the calculations. The Gomez et al. (2010) study was focused on WAIS collapse, and they highlighted the importance of incorporating shoreline migration in predictions of sea level fingerprints following the retreat (or advance) of grounded ice in marine settings. Finally, Mitrovica et al. (2011) summarized various technical issues associated with the calculation of sea level fingerprints, investigated the relative accuracy of published predictions, and demonstrated that 3D variability in mantle elastic structure would have negligible impact on these predictions.
Nearly all predictions of sea level fingerprints have, to date, assumed that the melting events are sufficiently rapid that the model for Earth’s mantle adopted in the calculations may be treated as being purely elastic. The exceptions are Bamber et al. (2009) and Mitrovica et al. (2009), who explored the importance of viscous effects in multicentury simulations of melting from WAIS and concluded that these effects were relatively minor. Both studies adopted one-dimensional (1D) mantle viscosity profiles common in global ice age studies, characterized by values in the range from 5 × 1021 to 1022 Pa s. These models do not capture the complexity of the mantle’s viscoelastic structure beneath WAIS, including both the range and lateral variability of viscosity.
The WAIS sits atop a major continental rift system (Wörner 1999) that is thought to have experienced at least two phases of extension since the breakup of Gondwana (e.g., Cande et al. 2000). Seismic evidence indicates that the lithosphere in the region is relatively thin (~65 km), in contrast with the stable East Antarctic craton, which has a thickness in excess of about 200 km (e.g., Morelli and Danesi 2004; Heeszel et al. 2013). In addition, anomalously slow seismic velocities characterize the sublithospheric mantle (to a depth of 200 km) below West Antarctica (Ritzwoller et al. 2001; Morelli and Danesi 2004; Hansen et al. 2014; Lloyd et al. 2015; Heeszel et al. 2016). Geodynamic modeling of the geological record in the area (Faccenna et al. 2008; Austermann et al. 2015) and seismic imaging (Emry et al. 2015) both indicate that this seismically slow region of the mantle is characterized by an active, thermal upwelling. Taken together, these studies suggest significant, 3D variations in lithospheric thickness and upper-mantle viscosity, with a particularly strong gradient as one moves across the shallow mantle beneath the western to the eastern sectors of the continent.
Recent models of ongoing glacial isostatic adjustment (GIA) across Antarctica in response to the last ice age have incorporated 3D variations in mantle structure (e.g., Kaufmann et al. 2005; A et al. 2013; van der Wal et al. 2015), with viscosities as low as 3 × 1018 Pa s beneath West Antarctica (Kaufmann et al. 2005). In this study we explore, for the first time, the impact of this variable viscoelastic structure on predictions of sea level fingerprints of relatively rapid WAIS mass loss. We will consider two case studies: mass loss over time scales of decades, characterized by thinning of ice cover, and centuries, where major collapse of grounded ice in marine-based settings is predicted (e.g., DeConto and Pollard 2016; Joughin et al. 2014; Levermann et al. 2014). The Maxwell time for mantle material with viscosity about 1019 Pa s is of order several years, and thus viscous effects may play a significant role in sea level changes over both time scales. That is, published fingerprints of WAIS melt based on elastic Earth modeling may be subject to significant error, particularly in the near field of the ice loss. The goal of this paper is to provide updated predictions of these fingerprints of melting in this region of complex, 3D viscoelastic Earth structure.
Our sea level projections are based on the gravitationally self-consistent sea level theory described by Gomez et al. (2010), which incorporates time-varying shoreline geometries and Earth rotation variations. The calculations are simplified relative to typical ice age sea level calculations because the initial topography is known at the outset. In applying the sea level theory, viscoelastic adjustments of the solid Earth in response to both surface mass (ice plus ocean) loading and perturbations in the centrifugal potential are computed using the finite-volume software described in Latychev et al. (2005). Load-induced polar motion is calculated following the theory of Mitrovica et al. (2005).
Recent updates to the finite-volume software allow regional grid refinement within a global model domain of lower resolution. In the calculations described below, the zone of grid refinement extends from the South Pole up to approximately 60°S (the outer boundary is not axisymmetric) and vertically down from the surface to 350-km depth; the grid is characterized by an average spatial (horizontal and vertical) resolution of 6 km to the base of the Moho, 12 km to a depth of 220 km, and 25 km to 350-km depth. Outside this region, these numbers increase to 12, 25, and 50 km, respectively, where the latter resolution extends to the base of the mantle.
The sea level calculations require two main inputs: a model for 3D viscoelastic mantle structure and the space–time history of ice cover. We describe each of the inputs, in turn, below.
a. Viscoelastic Earth models
All Earth models have elastic and density structure given by the 1D seismic Preliminary Reference Earth Model (PREM; Dziewonski and Anderson 1981) and assume a Maxwell viscoelastic mantle rheology that is compressible in the elastic limit. We will consider three different Earth models distinguished by their structure within the mantle. The first is a purely elastic model (MEL). In this case, the computed sea level change depends only on the net change in ice volume between the beginning and end of the calculation and not on the duration of the simulation. (In practice, the elastic sea level change is computed by melting all ice in a single, instantaneous time step.) As noted above, the vast majority of fingerprint calculations have been performed using purely elastic Earth models. The second Earth model (MVE1D) assumes a depth-varying viscosity profile comprised of an elastic (effectively infinite viscosity) lithosphere of 96-km thickness and uniform upper- and lower-mantle viscosities of 5 × 1020 and 5 × 1021 Pa s, respectively. The model falls within the class of 1D viscosity models favored in most GIA-based inferences of mantle viscosity (e.g., Lambeck et al. 1998; Mitrovica and Forte 2004).
Finally, a third Earth model (MVE3D) incorporates 3D variability in both the elastic thickness of the lithosphere and mantle viscosity. For the former, we adopt the lithospheric thickness model by An et al. (2015a) for the Antarctic plate and that by Conrad and Lithgow-Bertelloni (2006) for the remaining plates. Both models are scaled to produce a mean lithospheric thickness of 96 km, which results in a mean thickness of 65 km across West Antarctica and 200 km across East Antarctica (Fig. 1a).
The Earth model MVE3D is constrained to have a spherically averaged depth profile that matches the model MVE1D. The 3D viscosity variation below the lithosphere is estimated from seismic velocity heterogeneity using the multistep method discussed by Austermann et al. (2013). Because of the different geographical coverage and differences in resolution, we used a combination of three different seismic S-wave velocity models: within West and central Antarctica we used the model by Heeszel et al. (2016), for East Antarctica and the tip of the Antarctic Peninsula we used the model by An et al. (2015b), and for the rest of the globe we adopted the model called S40RTS (Ritsema et al. 2011). Figure S1 of the supplemental material illustrates these different domain geometries. The first two of the seismic models extend to 350-km depth, and S40RTS is adopted at greater depths across the entire mantle.
We derived a regional average seismic velocity model by moving sequentially through the three seismic models, beginning with the S40RTS model (Ritsema et al. 2011). In particular, the S40RTS model was adopted as the global model. The An et al. (2015b) model was then patched into this global model, and at each depth its velocities were shifted up or down so that the average value across the full Antarctic domain of the An et al. (2015b) model matched the average of the S40RTS model over the same Antarctic domain. The Heeszel et al. (2016) model was then patched into the An et al. (2015b) model so that at each depth the average across its smaller domain matched the average across the same domain of the An et al. (2015b) model. Variations in seismic wave speed are mapped into temperature perturbations following the depth-dependent scaling profiles described in Latychev et al. (2005). For the scaling of temperature anomalies to viscosity perturbations we used the scaling factors 0.037, 0.023, and 0.04 for the models by Heeszel et al. (2016), An et al. (2015b), and S40RTS, respectively [see Austermann et al. (2013) for the definition of this parameter]. The scaling factor varies in order to account for the amplitude differences between the S-wave models and thus to avoid significant viscosity discontinuities across the model domains. The resulting viscosity model is shown in Figs. 1b–f and results in viscosities as low as about 1018 Pa s below Marie Byrd Land (blue asterisk in Fig. 1a).
b. Ice melt geometries
We will consider two different WAIS melt geometries. The first (WA1), which follows most previous fingerprint analyses, will remove a thin uniform layer of ice across WAIS. While this is not a realistic melt scenario, our goal is to focus on the geographic perturbation to the sea level fingerprint driven by Earth model variability in the absence of changes in this signal associated with ice melt geometry. We will consider three different cases distinguished on the basis of the melt duration: 25, 50, and 100 yr.
The second, more extensive melt geometry (WA2) is taken from a coupled ice sheet–sea level simulation of future WAIS melting described in Gomez et al. (2015). Figure 2 shows snapshots of ice thickness at the start, middle, and end of the simulation (we do not include any melt of the East Antarctic Ice Sheet). The simulation had a duration of 1000 yr, but for our purposes we simply treat it as a realistic scenario for the geometry of WAIS collapse and treat the total duration of the collapse as a free parameter in the sea level fingerprint modeling. In particular, we additionally consider durations of 100, 300, and 500 yr.
All of our sea level predictions will be normalized by the total eustatic sea level change associated with the melt scenario and are therefore nondimensional. The eustatic sea level change is the rise in sea level computed for each melt scenario by assuming no deformation of the solid surface or self-gravitation of the surface mass load. For the WA1 melt scenario, the eustatic sea level change is 1 mm. In the case of scenario WA2, the effective eustatic sea level change (i.e., the uniform sea level change once all exposed marine-based sectors are filled with meltwater; Gomez et al. 2010) is 3.6 m, consistent with estimates of maximum ice mass loss from WAIS (Bamber et al. 2009). As discussed by Mitrovica et al. (2011), normalized sea level fingerprints can be scaled to consider melt events of similar geometry but different net changes in ice mass.
3. Results and discussion
Our goal is to quantify the potential signal from viscous effects in sea level fingerprints of rapid WAIS melting, and therefore results obtained for viscoelastic Earth models MVE1D and MVE3D will be shown relative to their associated (normalized) elastic fingerprint (Fig. 3). These elastic calculations show the features characteristic of previously published (elastic) sea level fingerprints of WAIS melting, including a zone of sea level fall in the vicinity of the ice melt that peaks about 9.4 and about 17.5 times the eustatic value for the WA1 and WA2 melt events, respectively, and a far-field sea level rise that peaks about 26% (for WA1) and 37% (for WA2) higher than the eustatic value in the North Atlantic, North Pacific, and Indian Oceans. The amplified signal in the near field is due to the more localized melt geometry. In the far field, the amplified signal for the WA2 simulation, relative to the WA1 case, is associated with an additional contribution in the former associated with meltwater being displaced by the uplift of exposed, marine-based sectors (Mitrovica et al. 2009; Gomez et al. 2010).
Figure 4 shows results in the Antarctic region for predictions based on the uniform ice melt geometry WA1 and three simulations distinguished on the basis of the duration of the melt event: 25, 50, or 100 yr. In all cases, we plot the total sea level change from the start to the end of the deglaciation. We note that if we continued the simulation beyond the end of the deglaciation, the impact of viscous effects on the sea level predictions would increase, and in this sense the results in Fig. 4 represent lower bounds on this impact.
We begin by considering viscoelastic calculations based on the 1D viscosity profile MVE1D (Fig. 4, left). As one would expect, the discrepancy from the elastic simulation (i.e., the signal from viscous effects) increases as the time scale of the melt event increases; specifically, the peak difference from the elastic case (Fig. 3a, inset) over West Antarctica is 1.1%, 2.1%, and 4.5% for the 25-, 50-, and 100-yr simulations, respectively. As we have noted, the Earth model MVE1D is within the class of 1D Earth models commonly inferred in GIA analyses. We conclude that the predicted sea level fingerprint of WAIS melt based on such models will exhibit only minor differences from the associated elastic fingerprint for time scales of ice melt on the order of a century.
Next, we consider analogous results based on the more realistic viscoelastic Earth model MVE3D (Fig. 4, right). In contrast to predictions based on model MVE1D, viscous effects become pronounced within just a few decades when MVE3D is adopted. In the near field, these viscous effects are characterized by localized sea level fall resulting from postglacial uplift in West Antarctica and an offshore zone of localized sea level rise resulting from the subsidence at the periphery, both of which are common features of ice age sea level predictions (e.g., Milne and Mitrovica 2008). The amplitude of these effects is nonnegligible; for the 25-, 50-, and 100-yr melt simulations, the peak sea level fall due to viscous effects occurs over the zone of minimum mantle viscosity (Marie Byrd Land; see Figs. 1b–f) and reaches amplitudes of 22.8%, 35.4%, and 54.1% of the peak elastic signal. Thus, predictions of the near-field sea level fingerprint of WAIS melting spanning the next few decades to century should be based on modeling that incorporates the complex mantle viscosity structure below the region.
In Fig. 5, we turn to the simulations of WAIS collapse (model WA2) and consider collapse duration time scales of 100, 300, 500, and 1000 yr. It is highly unlikely that a collapse of the magnitude shown in Fig. 2 could occur in the next century (the maximum rate of eustatic sea level change is ~50 mm yr−1 in this case), but this case serves as a bridge to the time scales considered in the uniform melt case of Fig. 4. We again note the increased amplitude of localized zones of postglacial uplift (sea level fall) and peripheral subsidence (sea level rise) driven by viscous flow as the adopted duration of WAIS collapse increases. Moreover, as in the results of Fig. 4, predictions based on the 3D model MVE3D significantly underestimate the impact of viscous effects on the WAIS fingerprint. In particular, while model MVE1D predicts a peak sea level fall over West Antarctica that is 5.6%, 16.5%, 26.9%, and 50.5% greater than the elastic fingerprint for melt durations of 100, 300, 500, and 1000 yr, respectively (Fig. 5, left), the analogous amplifications predicted using the 3D viscoelastic Earth model MVE3D are 53.1%, 110.5%, 152.3%, and 222.2% (Fig. 5, right).
To emphasize the importance of the viscous signal in the near field of WAIS, Fig. 6 shows the normalized elastic fingerprint and the normalized total viscoelastic fingerprint computed using the realistic Earth model MVE3D for two cases: ice melt scenario WA1 and a melt duration of 50 yr (Fig. 6a,b) and the collapse scenario WA2 and a melt duration of 500 yr (Fig. 6c,d). In both cases, viscous effects clearly imprint a significant signature on the pattern and magnitude of the computed sea level fingerprints throughout the near field.
As a final numerical test, we considered the sensitivity of the results to variations in the magnitude of the lateral variability in mantle viscosity. As we noted above, the minimum viscosity below West Antarctica in the case of Earth model MVE3D is approximately 1018 Pa s. We ran two additional simulations in which the scaling from seismic velocity to viscosity was varied to yield minimum viscosities of 2 × 1017 and 5 × 1018 Pa s. These calculations were based on the WAIS collapse scenario WA2 and a melt duration of 500 yr. Figure 7 shows the total viscous signal (i.e., full calculation minus elastic results) for these simulations, together with the analogous result for model MVE3D (reproduced from Fig. 5, bottom right). As one would expect, the viscous component of the sea level response increases with decreasing mantle viscosity; however, this signal remains significant even for the prediction based on the highest-viscosity Earth model. In particular, in order of decreasing viscosity, these simulations yield a maximum sea level fall that is 109.4%, 152.3%, and 193.6% of the elastic response, respectively, at the end of the melting event.
Fingerprints of sea level change following the rapid melting of an ice sheet are almost exclusively computed by assuming that the solid Earth response to melting is purely elastic. This assumption is questionable in the case of ice mass changes in West Antarctica given that geological and geophysical evidence suggests that sublithospheric viscosity below the region is approximately two orders of magnitude lower than values typically adopted in ice age sea level calculations. We have demonstrated, using a 3D viscoelastic Earth model consistent with available constraints, that sea level predictions in the near field of WAIS will depart significantly from previously published elastic fingerprints in as little as a few decades. This result has implications for studies of sea level change due to both ongoing mass loss from WAIS over the next century and large-scale collapse of WAIS on centennial-to-millennial time scales.
We have also demonstrated that the amplitude of these viscous effects is not accurately modeled using standard 1D viscosity profiles. As an example, both Bamber et al. (2009) and Mitrovica et al. (2009) considered the impact of viscous effects on the elastic fingerprints of WAIS collapse over time scales of centuries and longer by adopting 1D Earth models very similar to MVE1D. We conclude that their calculations significantly underestimated the impact of viscous relaxation in the near field of the sea level fingerprints. More generally, our results also suggest that inferences of Antarctic ice volumes since the Last Glacial Maximum based, in part, on analyses of sea level datasets from the region and models akin to MVE1D (e.g., Whitehouse et al. 2012) may need to be reevaluated.
Finally, the strong and geographically variable near-field signal of viscous effects associated with the 3D model MVE3D (Fig. 6) suggests that coupled ice sheet–sea level models of WAIS evolution (e.g., Gomez et al. 2013, 2015) should, in future work, incorporate this complexity in viscoelastic Earth structure.
This work was supported by the U.S. National Science Foundation (Grants ARC-1203415, OCE-0825293, PLR-1142518, and PLR-1246712), the NJ Sea Grant Project 6410-0012, the Natural Sciences and Engineering Research Council of Canada Discovery Grant program, the Canada Research Chairs program, Harvard University, and McGill University.
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-16-0388.s1.