## Abstract

The Antarctic mantle and lithosphere are known to have large lateral contrasts in seismic velocity and tectonic history. These contrasts suggest differences in the response timescale of mantle flow across the continent, similar to those documented between the northeastern and southwestern upper mantle of North America. Glacial isostatic adjustment and geodynamical modelling rely on independent estimates of lateral variability in effective viscosity. Recent improvements in imaging techniques and the distribution of seismic stations now allow resolution of both lateral and vertical variability of seismic velocity, making detailed inferences about lateral viscosity variations possible. Geodetic and palaeosea-level investigations of Antarctica provide quantitative ways of independently assessing the 3D mantle viscosity structure. While observational and causal connections between inferred lateral viscosity variability and seismic velocity changes are qualitatively reconciled, significant improvements in the quantitative relations between effective viscosity anomalies and those imaged by P- and S-wave tomography have remained elusive. Here we describe several methods for estimating effective viscosity from S-wave velocity. We then present and compare maps of the viscosity variability beneath Antarctica based on the recent S-wave velocity model ANT-20 using three different approaches.

Solid Earth rheology is that discipline which deals with the laws that govern the deformation of rocks under the high-temperature and -pressure conditions of the mantle. Observations linked to rheology require input from many disciplines in the geological and geophysical sciences, and examples include composition, temperature and pressure conditions from xenoliths, elastic parameters from seismic wave propagation, and density variations derived from gravity, topography and magnetic anomaly data (e.g. Scheinert *et al.* 2016; Martos *et al.* 2017; Pappa *et al.* 2019). Simulations of mantle circulation (Bredow *et al.* 2021, this volume), glacial isostatic adjustment (GIA), post-seismic deformation (Barletta *et al.* 2022, this volume) and their surface expressions require quantifying large variations in mantle rheology. For several areas of the world we have some quantitative bounds on lateral variability in upper mantle viscosity, such as the contrast between Fennoscandia and Iceland (Wang *et al.* 2008). Antarctica has long been known to have contrasting mantle and lithosphere properties between east and west, somewhat akin to the contrasts between the Precambrian craton of eastern Canada and the geologically younger crust and mantle of western North America (Crittenden 1963; Austermann *et al.* 2020).

A broader goal of this chapter is to make quantitative connections between flow laws derived in the rock deformation laboratory and the integrative power of the datasets for constraining crustal motion and Earth structure using a network of global positioning system (GPS) and broadband seismic instruments. A period of intense data collection began after the 4th International Polar Year of 2007 and was organized as the POLENET project. Much of this integrated dataset may be brought to bear on understanding the nature of the Antarctic lithosphere (e.g. Nield *et al.* 2018; Shen *et al.* 2018; Scheinert *et al.* 2021, this volume). An overview of seismic imaging of the Antarctic crust, lithosphere and upper mantle is given in the chapter by Wiens *et al.* (2021, this volume) that includes tectonophysical interpretations, thus complementing the study of Earth rheology. The reader may also find the recent review by Jordan *et al.* (2020) to be useful. However, in this chapter we focus on the mantle response to changes in surface ice loading over centennial to millennial timescales with steady-state high-temperature creep. The rate of strain associated with this creep may be causally connected to both vertical uplift and horizontal crustal motion rates that have been measured using bedrock GPS (e.g. James and Ivins 1995; Groh *et al.* 2012; Barletta *et al.* 2018). This motion also occurs where it cannot be measured, deep below the surface of the ice sheet, over extensive regions of the vast ice-covered interior. Such extensive continental-scale mantle and crustal motion is important to model since the vertical component creates a time-varying gravity field. The modelled GIA signal needs to be removed from gravimetry measurements, such that accurate ice mass-balance trends for the Antarctic ice sheet and its contribution to ongoing sea-level rise may be retrieved using space gravity mission data (e.g. Caron and Ivins 2020).

Early surface-wave studies in Antarctica (e.g. Knopoff and Vane 1978; Danesi and Morelli 2001; Ritzwoller *et al.* 2001) revealed the existence of a major contrast in the structure of the uppermost 300 km of the mantle between East and West Antarctica. The boundary between these deep structural elements is roughly coincident with the surface expression of the Transantarctic Mountains (TAM: see Fig. 1). The seismic data retrieved since 2007 are of sufficiently high quality that we can begin to quantify the lateral heterogeneity in the viscous response of the mantle to ice-mass surface-load changes. Lateral heterogeneity in the mantle is important for reconciling observed and predicted crustal motion in GIA models (Powell *et al.* 2020), and is an important aspect for determining the stability of the ice sheet over a changing bed topography, both in the past and in the future (Gomez *et al.* 2018). Lateral variability in elastic-wave velocity (see Fig. 1) is related to lateral temperature variations. Using the temperature-dependence of elastic constants and anelasticity determined by theory and experiment, we can try to relate seismic velocity change to changes in effective mantle viscosity, since the latter are known to be strongly dependent on temperature.

The focus of this chapter is to examine the feasibility of developing causal interconnections between seismic tomographic imaging, mantle viscosity retrieved from Maxwellian 1D GIA models and flow laws for steady-state creep quantified by high-pressure and -temperature laboratory experiments. This is a rather grand challenge, and may have consequences for regions other than Antarctica. However, here we restrict our study only to a single tomographic S-wave imaging study (Lloyd *et al.* 2020) for Antarctica and restrict ourselves to the upper mantle. The feasibility of establishing causal interconnection is tested by conducting an intercomparison of three methods at three different depths. The methods for each are documented in the remaining pages of this chapter. In Table 1 we give a short synopsis of some of the differences in the methods. Later we refer to these as Approaches 1, 2 and 3 (see the first column in Table 1). Not only is the selection of laboratory flow laws generally different, so are the implementation steps yielding anelastic slowing of the shear-wave velocity. The latter are necessary to derive the appropriate thermal deviation *δT*(*r*, *θ*, *ϕ*) from tomography. We anticipate from the outset that there will be large differences between the predictions. These should manifest themselves as differences in mapped log_{10}Δ*η*(*r*, *θ*, *ϕ*) for the Antarctic mantle and in comparisons to the inferences from 1D GIA models compatible with geodetic data.

The scope of study in this chapter is limited by the fact that Antarctica is an ice-sheet-covered continent with a logistically difficult environment in which to obtain bedrock GPS, relative sea-level (RSL) data and maintain seismic station operation. We are perhaps blessed by the fact that geochronological and ice-core constraints on the ice-sheet history during the late-Quaternary place tight constraints on the millennial timescale ice-load history (e.g. Bentley *et al.* 2014; Anderson *et al.* 2014; Whitehouse *et al.* 2017), a critical ingredient for establishing viable GIA predictions. Neither a comprehensive 3D compressional (P-wave) map can be combined with S-wave mapping (e.g. Goes *et al.* 2000), nor does an attenuation map of the upper mantle for Antarctica exist (e.g. Lawrence and Prieto 2011). This fact diminishes the capability of retrieving higher-fidelity mapping of *δT*(*r*, *θ*, *ϕ*) that can be produced in upper mantle environments where seismic station coverage is denser in both space and time (e.g. Goes *et al.* 2000). Nonetheless, in this chapter we demonstrate that a coherent rationale exists for quantitatively connecting Antarctic seismic tomography to GIA-determined mantle viscosity and laboratory flow laws. This inference may have an impact on the construction of future models of glacial isostatic adjustment with lateral variation in Earth structure (GIA-LV) and for formulating paradigms for their uncertainty.

In this chapter, we will first review the factors controlling mantle rheology as well as seismic-wave velocities in the mantle. We then discuss a general computational strategy for estimating 3D mantle viscosity from seismic observations and present new maps of viscosity for the Antarctic upper mantle derived in three different ways. In our conclusions, we attempt to draw out robust predictions that are common to the three methods. A common starting point for our estimates is the ANT-20 S-wave velocity model (Lloyd *et al.* 2020). This model takes advantage of recent increases in seismic station density and computational techniques to derive higher-resolution images of mantle structure using full waveform adjoint tomography (Wiens *et al.* 2021, this volume). The focus of this chapter is on the most basic elements of the large uncertainties encountered in deriving viscosity from seismic models. We identify and detail a number of the relevant quantitative elastic and plastic deformation parameters, but assess their relative uncertainty in a more qualitative way. Our main conclusions will, hopefully, offer clarity and guidance for further efforts to reduce this uncertainty, and thereby improve our estimation of mantle rheology.

## Background

Deformation in the mantle is controlled by the mobility of defects in the crystal structure, either within or at the boundaries of rock grains (Kohlstedt and Hansen 2015; Masuti *et al.* 2019). Therefore, this chapter begins with a brief synopsis of some of the microphysical models derived from laboratory deformation experiments. Such experiments reveal the same defect structures that are found in field samples of mantle rocks (e.g. Matsyuk *et al.* 1998; Johanesen and Platt 2015).

While the essence of rock flow at mantle temperature and pressure is to be found microscopically, all geodynamic simulations rely on the continuum hypothesis. This hypothesis asserts that the constitutive law governing the relations between stress and strain, and their time derivatives, are independent of scale. While a computational convenience is to describe mantle flow controlled by a constitutive equation that has an effective Newtonian viscosity, a large body of experimental results suggests a non-Newtonian behaviour. Hence, any discussion of Antarctic mantle rheology must probe this issue. Detailed descriptions of the experimental basis of the laws derived from studies of the mechanical strength of mantle rocks at high temperature and pressure is well beyond the scope of this chapter.

Overall, there exists a gap between laboratory-based flow laws and the rheology obtained by models of geodetically measured solid Earth deformation. This gap may only be reconciled by placing both methodologies into a rigorous statistical framework that will define a subset of the flow laws and model predictions that are mutually compatible (e.g. Jain *et al.* 2019). It is hoped that by writing this chapter we can lay some groundwork that might further this aim for the Antarctic mantle.

### The role of experiments

Figure 2 shows a classic representation of the deformation diagram used in study of mantle rheology to define regimes of control by microscopic mechanisms as they vary with stress, strain rate, temperature and/or grain size. Here we show one example of a deformation map produced by Karato (2010) for olivine, a rock type that dominates the mantle above 410 km depth. The portion of the deformation diagram in Figure 2 that shows stress values below 100 MPa corresponds to the regime in which GIA-generated stresses occur. GIA stresses may reach 10 MPa in Antarctica near the lithosphere–asthenosphere boundary for realistic ice loading and unloading (e.g. Ivins *et al.* 2003). In this regime, two microscopic mechanisms may control effective rheology: diffusion or dislocation creep. Later in this chapter, both of these mechanisms are to play a role in the scaling to S-wave seismic tomography.

There are caveats that are well known when applying experimentally derived flow laws to geodynamic and GIA models. An obvious caveat is the difference between mantle strain rates (𝒪 10^{−15} − 10^{−12} s^{−1}, where 𝒪 indicates the ‘order of magnitude’) and those of the laboratory (𝒪 10^{−6} − 10^{−3} s^{−1}) (see the light contour lines in Fig. 2). There is the additional caveat that the principal shear stresses are higher in the laboratory (see Fig. 2). Some experiments at higher strain rate (𝒪 10^{−3} s^{−1}) and higher stress (𝒪 80 − 100MPa) are relevant to shear localization in the lithosphere (e.g. Hansen *et al.* 2012). In Figure 2, the upper portion shows a region of the deformation map (deviatoric stress v. rock grain size) that corresponds to the material response that involves a plastic yield stress. This type of deformation may occur at relatively low temperature and high stress. In this regime, on a microscopic level, dislocation glide may dominate, a mechanism generally known as Peierls mechanism. This mechanism may be important to deep lithospheric shear zones (Kameyama *et al.* 1999). Such rheology has not been accounted for in GIA models and we do not treat this rheology in this chapter.

Constraints provided by seismology can be linked to temperature- and pressure-dependent creep laws for mantle rock aggregates composed of typical mantle minerals (olivine, pyroxene and Al-silicate). A modicum of understanding of the physics involved in both elastic wave propagation *and* dislocation and impurity diffusion within the rock grain interior and at grain boundaries is required. These will be briefly reviewed in this paper, along with some of the caveats that accompany these linkages. Furthermore, we assemble the current state-of-the-art interconnections between the high-pressure and -temperature thermodynamics and solid-state relations that allow us to quantify the effective mantle viscosity.

### Steady-state v. transient laboratory experiments

The experimental set-ups for steady-state and transient creep are described by Paterson and Olgaard (2000), Karato (2010) and Kohlstedt and Hansen (2015). Two fundamental quantities must be extracted from such experiments for our application. These are the activation enthalpy for steady-state and transient, or anelastic, creep, *H** and *et al.* 1982; Jackson 2019). Attaining a true steady state in the static experiments, wherein the rate of creation and destruction of microscopic elements at grain boundaries and within grains reaches equilibrium, is one of the challenges for developing creep laws that are acceptable in modelling mantle flow. As a rock deformation experiment begins, those microscopic elements that are capable of being mobilized under an imposed shear stress move relatively quickly. After a time, these elements (e.g. gliding or climbing dislocations and grain-boundary slip) may become blocked by obstacles, including their own bending or tangling, and, as a consequence, the creep rate of the rock aggregate slows. Technically, the macroscopic phenomenon is referred to as strain hardening (e.g. Peltier *et al.* 1980) and formal connections may be made between the damping of seismic waves to the microscopic viscoelastic strain (e.g. Minster and Anderson 1980; Jackson 2019). The mechanisms controlling the transition from transient to steady-state creep are numerous (e.g. Hanson and Spetzler 1994). While complex, these are linked to the anelastic properties of mantle rocks (e.g. Faul and Jackson 2015; Jackson 2019; Masuti *et al.* 2019).

### Some general geophysical observations

Although our understanding of the physics of creep at upper and lower mantle conditions is limited, there is a clear observational connection between regions of low seismic-wave velocity and rapid isostatic response measured by geodetic uplift and/or rapid post-glacial land emergence (e.g. Sigmundsson 1991). While accounting for lateral heterogeneity in mantle deformation has been addressed over many decades (Wu 1999), today both higher-resolution seismic tomography and more dense and accurate geodetic networks are reporting data, meaning that the study of lateral heterogeneity in mantle rheology is evolving into a new frontier in solid Earth sciences (e.g. Austermann *et al.* 2013; Lau *et al.* 2018; Huang *et al.* 2019). Estimates of viscosity based on seismic-wave speeds and high-temperature creep experiments can be tested using geodetic data that are capable of constraining GIA or post-seismic modelled crustal motions. This enhances our knowledge of regional mantle rheological structure, and produces more confident viscosity estimates than in continental areas that lack such geodetic data constraints.

### Scaling S-wave velocity to viscosity: retrospective for GIA

The general study of GIA has a substantial impact on how we interpret the rheology of the Earth's mantle and lithosphere over millennial to million year timescales. The preponderance of the GIA models that are used to simulate gravitational viscoelastic flow over the ice ages assume a radially stratified mantle that is approximated as a 3D linear Maxwell viscoelastic material in tensor form (e.g. Peltier *et al.* 1981). In this model rheological framework, the Earth responds viscously and elastically over long and short timescales, respectively. Its success in giving rational explanation to both geodetic (e.g. Milne *et al.* 2004) and RSL data (e.g. Lambeck *et al.* 2014) has meant that extensions to Maxwellian set-ups that contain non-linear power-law viscosity (e.g. Gasperini *et al.* 2004) generally fall outside of the GIA model paradigm most commonly used in cryospheric sciences, by geodesists and in the palaeo-oceanography community (e.g. Wahr *et al.* 2000; Whitehouse 2018; Dobslaw *et al.* 2020).

The structural roots of modelling GIA with laterally varying Earth structure (GIA-LV) have a number of commonalities. Seismic tomography provides elastic-wave velocity in 3D space (*r*, *θ*, *ϕ*), where the three spatial variables are radial position, co-latitude and longitude, respectively. A systematic way to describe these data is in the form *δν*_{S}, represents the deviation of velocity from a mean value, *δν*_{S} are related to absolute temperature. The temperature-dependence of the elastic constants governing wave propagation are known through theoretical and semi-empirical methods (e.g. Kumazawa and Anderson 1969). All of the GIA-LV models developed to date, including those that incorporate non-linear power-law viscosity, have assumed that the deviations from the mean are caused by a temperature deviation *δT*(*r*, *θ*, *ϕ*) away from a radially dependent mean *T*_{0}(*r*). Paulson *et al.* (2005), for example, assumed that the deviations in *δν*_{S} were linearly related to lateral density anomalies *δρ*, such that:
*et al.* (2005) also assumed that the relationship of seismic velocity to temperature is completed by additionally assuming that equation (1) is divided by the thermal expansivity, *α*_{th}; thus, deriving the spatial dependence of the thermal deviation *δT*(*r*, *θ*, *ϕ*) from the mean *T*_{0}(*r*). The reader may note that this is the method adopted by A *et al.* (2013) for GIA-LV predictions for Antarctica. Both a nominal radial average temperature and its aspherical deviation were assumed to be known, and are written in the form:
*A*_{0}, a dimensionless activation energy, *g*, and *T*_{m}(*r*) is the melting temperature (e.g. Poirier 2000). Thus, all of the ingredients for performing finite-element model (FEM) solutions of the time-dependent sixth-order partial differential equation system for GIA-LV are defined.

Latychev *et al.* (2005) also developed a numerical finite-volume method for GIA-LV that also used equations (1) and (2), but made no assumption about the details of the governing rheology other than the Arrhenius-dependency given by the temperature in the exponential term. Therein, a viscosity ratio is defined as:
*ɛ* is used to adjust the total lateral variability in viscosity assumed in any GIA-LV computer simulation. A typical value for this adjustment parameter is *ε* ≈ 0.4 K^{−1} (e.g. Austermann *et al.* 2013). Like Paulson *et al.* (2005), Latychev *et al.* (2005) assumed:

Latychev *et al.* (2005) and Paulson *et al.* (2005) developed simple and flexible scaling relations, but which are not closely tied to rheological laws that are currently in debate in the experimental rock deformation community (e.g. Kohlstedt and Hansen 2015). In contrast, the FEM formulation for GIA-LV by Wu (2005) sought to expand the scaling method so that it explicitly accounted for activation enthalpy of high-temperature creep, its decomposition into activation energy and activation volume, along with the pressure dependence. This was accomplished by adopting the formulation of Ivins and Sammis (1995). The latter method also forced the globally averaged radially-dependent viscosity, *η*_{0}(*r*), to be defined by the same rheological dependencies. This guarantees an internal self-consistency, such that the rheological parameters could be tested against the most robust among competing 1-D GIA inversions for mantle viscosity. The formulation proposed by Ivins and Sammis (1995) used the same assumptions for the relationship of elastic shear-wave velocity to density and temperature as in equations (1), (2) and (5), although it more rigorously accounted for the Grüneisen parameter, a fundamental property of theoretical thermal elasticity (e.g. Isaak *et al.* 1992; Stacey and Hodgkinson 2019). These relationships will be explicitly given in the ‘Methods for Approach 1’ and ‘Methods for Approach 2’ sections later in this chapter.

Wu (2005) adopted the log_{10} representation of the viscosity field and employed the relationship to the activation enthalpy parameter, *H**, as proposed by Ivins and Sammis (1995). This representation can be written as the ratio:
*a* ≡ *H**/*R*_{G} and *R*_{G} being the universal gas constant. Upon taking log_{10} of both sides of equation (6), this becomes:
*δT* = 0:

By truncating the expansion at a single term, we arrive at the form employed in Wu (2005) and to be used in our intercomparison of methods. The approximation is also inherent to the strategy of Latychev *et al.* (2005). In Appendix A of this chapter, we use variability up to the maximum temperature differences found in the Antarctic thermal model of An *et al.* (2015) at a depth of 240 km to give a quantitative evaluation of this approximation. Although Wu (2005) used this form, a later GIA-LV model (Wang *et al.* 2008) extended the relation for connecting *δν*_{S} and *δT* to include anelastic corrections that were earlier advocated by Karato (1993). It is this anelastically corrected form that is the starting point of our intercomparative study entertained later in this chapter.

## Microphysics

It has long been understood that deformation mechanisms in ceramic materials at high temperature are controlled by the diffusion of ionic impurities (point defects) or by slip due to dislocations (linear or planar defects) in the crystal lattice. Ceramic materials are distinguished by their covalent and ionic atomic bonding, usually involving one or more oxygen atoms. This intracrystalline bonding allows experiments to be performed on certain analogue materials with crystal structure such as rock materials at mantle temperatures and pressures, thus relevant to both thermal elasticity and creep. This section discusses the main mechanisms, as well as a generalized flow law and the main rationale used to approximate an ‘effective’ viscosity.

### Steady-state flow law and temperature-dependent anelasticity

The diffusion of point defects and the slip of dislocations are the microscopic manifestations of diffusion and dislocation creep, respectively. Both creep mechanisms are controlled at some stage by diffusion and, hence, by temperature. The strain rate, *σ*, is then exponentially dependent on the local mantle temperature, *T*, typical of all thermal activation processes (e.g. Kittel 2004; Cressler and Moen 2012). The form is:
*H** is the activation enthalpy associated with diffusivity in solids (Poirier 2000; Karato 2008), *R*_{G} the universal gas constant (8.314 × 10^{−3} kJ mol^{−1} K^{−1}), *n* is a semi-empirical exponent (1 ≤ *n* ≤ 5) and ^{−1}Pa^{−n}, which is often tagged to a reference strain rate *H** at a series of temperature values. In this chapter we shall make use of two activation enthalpies: one for steady-state creep and the other for anelastic and transient (viscoelastic) strain processes in mantle rock. These activation enthalpies are decomposed into two terms:
*P* is pressure, and pairs *E**, *V**, *E** is more straightforward to recover from laboratory experiments than values of *V**, and similarly for their transient counterparts, *V** must accommodate stepwise increases in pressure while carefully accounting for all other factors influencing creep mobility (Dixon and Durham 2018).

It is to be understood that these activation parameters are arguments in exponential functions multiplying *T*^{−1}. Hence, they represent the fact that there is a strong temperature-dependence both to steady-state creep and to seismic dissipation and dispersion (e.g. Wiens *et al.* 2008). The latter viscoelastic properties are essential for properly developing transfer functions between S-wave velocity and mantle viscosity.

The normal procedure for determining

If the constitutive laws only involved these two thermodynamic variables, while probably non-unique, the problem of mapping seismic velocity to viscosity using the methods of either Ivins and Sammis (1995) or Wang *et al.* (2008) is quite tractable, provided that both the temperature and the rate-controlling creep mechanism are decipherable. However, not only is it difficult to decipher the rate-controlling mechanism, but the pre-exponential constant *f* is the fugacity of the respective impurity in the rock material (Karato 2008). It also represents the partial pressure of the invasive atomic species in the crystal structure (here *T*(*r*, *θ*, *ϕ*). Therefore, they are briefly discussed in the following subsection.

### Role of water and grain size

Small grain sizes tend to promote material transport by grain-boundary diffusion (Karato and Wu 1993), whereas large grain size promotes deformation by dislocation creep. Also, experiments with mantle olivine and synthetic anorthosite indicate that trace amounts of structurally bound water have a pronounced weakening effect (Mackwell *et al.* 1985; Hirth and Kohlstedt 1996; Dimanov *et al.* 1999; Mei and Kohlstedt 2000), thus indicating the importance of water-weakening in increasing the mobility of rock in the lithosphere and upper mantle environment. A wet sample of rock is defined when there are roughly 500–3000 H per 10^{6} Si and a dry sample at one order of magnitude smaller fraction H content (1 ppm H/Si ≈ 10^{−5} wt% of water in olivine, or 5 × 10^{−3}–3 × 10^{−2} wt% H_{2}O). Proper quantification of the effect of water-weakening at mantle temperatures and pressures remains controversial (Fei *et al.* 2013). Generally, volatile impurities, such as water, are a subclass of point defects. Other atomic or molecular species may also be important in controlling mobility.

In order to account for grain size in the constitutive equations, *Ad*^{−p} (Karato *et al.* 1986; Dimanov *et al.* 1999), and for both grain size and wetness content to *A* are dependent on the inclusion, or exclusion, of the effects. Here *p* and *q* are experimentally determined exponents. Possibly the most complete prefactor would be to set:
*et al.* 2017) and water content (O'Donnell *et al.* 2017) might be important for Antarctica. Examining these latter causes is warranted should we discover especially poor resemblance of effective viscosity scaled to shear-wave velocity maps, *ν*_{S}(*r*, *θ*, *ϕ*).

It is worth noting that both structurally bound water and grain size are likely to affect seismic velocities, although quantitative relationships are uncertain and, thus, difficult to implement formally in viscosity estimates. Water is proposed to reduce seismic velocities by increasing anelasticity, through a similar mechanism to its effect on viscosity (Karato 2003). Recent laboratory experiments indicate a small direct seismic velocity perturbation due to the presence of water (Cline *et al.* 2018), yet a substantial sensitivity to oxygen fugacity,

### Deformation mechanisms in mantle conditions

Simulation of upper mantle conditions in the laboratory are limited by confining pressure, limited temperature ranges and relatively high strain rates that must be studied in the laboratory. Nonetheless, there are advantages in that experiments have been performed on natural or synthetic rocks of either the major constituents of the mantle or their analogues. Sampling of upper mantle rocks under a microprobe reveal those subgrain-scale dislocation sites that are associated with ceramic high-temperature creep mobility. If diffusion creep dominates deformation in the upper mantle, then *n* ≈ 1 in equation (9), an assumption often employed in GIA models. However, this requires the upper mantle to have small grain sizes throughout, contrary to observations (Karato and Wu 1993; Dannberg *et al.* 2017; Coltorti *et al.* 2021, this volume). For large grain size, or large stress, dislocation creep dominates, as in Figures 2 and 3. An alternative mechanism, which has been studied extensively over the past several years is grain-boundary sliding (GBS) (Hirth and Kohlstedt 2004; Hansen *et al.* 2011). While this mechanism may be controlled by diffusion for small grain size, experiments with naturally occurring olivine reveal that dislocation-accommodated GBS can dominate creep with *n* ≈ 2.9 at upper mantle conditions and realistic grain size (*d* ∼ 1–10 mm). Figure 3 shows the dominant deformation mechanisms in a stress–temperature diagram for cases of small (0.01 mm) and large (1 mm) grain sizes. As an example, consider the contour line corresponding to a strain rate of 10^{−15} s^{−1} for an olivine with a 1 mm grain size. At high temperature and low stress, diffusion creep will dominate. With increased stress, dislocation creep will become more important. The stress at which the transition occurs will be higher when grain size is decreased. For very small grain sizes (10 µm), diffusion creep will give way to GBS when stress is increased.

## Continuum approach and effective viscosity

In the continuum point of view taken in geodynamic flow simulations, material behaviour can be well described by the Maxwell model (e.g. Ivins *et al.* 1982; Karato 2008). In the Maxwell model, elastic deformation occurs instantaneously, while at relatively long timescales movement (≥1–10 kyr) is described by slow flow in a viscous fluid. GIA codes book-keep the elastic deformation, as this is important for RSL and Little Ice Age geodetic computation (Ivins *et al.* 2000; Sabadini *et al.* 2016), but they are usually ignored in geodynamic models wherein the total viscous strain vastly exceeds that of the elastic components over tens of millions of years, or greater. The Maxwell model is generally faithful to the microphysical processes; it corresponds to a situation in which defects can move indefinitely after an initial elastic strain (Karato and Spetzler 1990; Karato 2010). To use the Maxwell model, a viscosity parameter is needed. The viscosity in such models can be linked to both the laboratory-based laws and to seismically based mantle models by defining a so-called effective viscosity.

### Choice of effective viscosity relation

The advantage of defining an effective viscosity, *et al.* 2012) and to the boundary value problems typically employed in the GIA problem (Spada *et al.* 2011); three different forms that have been proposed by Poirier (2000) and Ranalli (2001). Each, however, require some flow variable to be fixed. The three relationships yielding an effective viscosity are given below. At constant local background strain rate:

Finally, at constant rate of energy dissipated in the flow:
*n* (e.g. see equation 12). The reader may consult Ranalli (2001) or Karato (2008, equations 19.3a–c) for the corresponding details. The Maxwell assumption so often used in GIA modelling is that either *η*, defined with the equation
*et al.* 2005).

Since the time of Norman Haskell's work in the 1930s, the simple Newtonian viscosity, *η*, is the sought after parameter for the mantle interior in various forward-inverse GIA and mantle-flow models. In what follows in this chapter, we will be using the effective viscosity as defined by equation (14). This will avoid employing the power-law exponent *n* for dislocation creep into the argument of the exponential in the flow law. We remark that in doing so, we choose the approximation that will tend to maximize the lateral variability predicted from the ANT-20 seismic tomography model.

### Maxwell models: why are they relatively successful?

Considering that most laboratory experiments yield *n* = 1) and dislocation (

## Viscosity estimates for Antarctica

As is clear from the previous section, viscosity involves many parameters that are themselves uncertain. Past studies of laterally varying viscosity beneath Antarctica (e.g. Kaufmann *et al.* 2005; A *et al.* 2013; van der Wal *et al.* 2015; Gomez *et al.* 2018) have all agreed on a larger East Antarctic viscosity, owing to the faster seismic velocities imaged there. The results of these studies, however, substantially differ in the predicted magnitudes. Although, the resolution of the mantle seismic structure has improved through the years, many of the differences may be caused by methods employed in transforming seismic velocity to viscosity.

In the several subsections that follow, we estimate Antarctic viscosity in three different ways. In the first approach, viscosity perturbations are calculated directly from seismic models according to the approach of Wu *et al.* (2013). These use the temperature derivatives of shear velocity from Karato (2008) to estimate a temperature anomaly and to estimate a viscosity anomaly relative to a reference viscosity model. The second method uses the structural form derived by Ivins and Sammis (1995), but is extended to capture the anelastic correction using the theoretical development of Karato and Karki (2001) and Karato (2008). The third method has similarities to the first two approaches, but assumes a complete representation of the rock experimental results, including the non-thermodynamic prefactors. Both the first and third approaches will rely on temperature derivatives as explicitly computed by Karato (2008). In all approaches, we assume that global average 1D seismic velocity and mantle temperature are models that are relatively well known. In addition, Approaches 1 and 2 assume that radial 1D viscosity models, *η*_{0}(*r*), for the upper mantle are well determined. They then compute a 3D viscosity model based on perturbations to these global averages based on a straightforward conversion from seismic velocity to effective viscosity. Naturally, in each case of Approach 1 and Approach 2, *η*_{0}(*r*) must also be logically linked to the flow law assumed. To give clarity to the formulations of the three methods, we give some of the requisite fundamentals of the classical thermodynamics and elasticity theory. To better elucidate the comparison and the uncertainties, the three methods all start from the same seismic structure model of Lloyd *et al.* (2020) (Wiens *et al.* 2021, this volume) termed ANT-20, which is the highest-resolution continental-scale model available. Figure 4 shows the ANT-20 *δ* ln*ν*_{S}(*r*, *θ*, *ϕ*) field at depths of 150 and 550 km. Together with the map in Figure 1, these show the lateral variations in shear-wave velocity common to all three approaches for estimating effective viscosity that will be discussed in the next four subsections of this chapter.

### Relations between elastic and seismic parameters

Seismic velocity in isotropic rock media is related to density, *ρ*, elastic bulk modulus, *κ*, and shear modulus, *μ*, by:

Thus, both shear (*v*_{S}) and compressional (*v*_{P}) seismic velocities are related to density. It is often useful to consider seismic velocities as anomalies relative to a 1D radial background (reference) model. For example, the change in shear velocity, *δν*_{S} (

In the following we focus on the shear-wave velocities, as shear-velocity models generally show the highest continent-wide resolution in the Antarctic upper mantle. Anderson (1987) derived relations from thermodynamic and solid-state physics principles that interrelate seismic velocity anomalies and density anomalies. One goal of linking global seismic mapping and geodynamical modelling is to constrain the ratio of variable seismic velocity to variable density, *δ*ln*ν*_{S}/*δ*ln*ρ* (e.g. Dziewonski *et al.* 1993). In theory, inversion of such ratios can be linked to thermodynamic quantities derived from pressure and temperature derivatives of elastic parameters of mantle mineralogy. For interpreting seismic velocity anomalies and converting them to viscosity, a fundamental interrelation used by Ivins and Sammis (1995) was:

The ratio appearing on the right-hand side of equation (21) is the Grüneisen parameter for the shear modulus, *c*_{1}:
*c*_{1} is identical to *δ*_{μ} as defined in the monograph by Karato (2008, equation 4.35b). The thermal Grüneisen parameter is defined as a ratio of vibrational frequencies to molecular volume, *γ*_{μ}:

The relations we have presented in this subsection establish the basic theory for relating seismic-wave velocity to the Grüneisen parameter, a fundamental property of the thermodynamic state of crystalline rock at high temperatures and pressures. In the next subsection we show how this is a key link to temperature derivatives of the shear-wave velocity. Here we strictly employ the thermodynamic and elasticity approach that is advocated by D. L. Anderson (1988), O. L. Anderson and Isaak (1995), Karato (2008), Stixrude and Jeanloz (2015), Stacey and Hodgkinson (2019) and others. We must also be cognizant of the uncertainties associated with the thermal elastic constants of mineral physics at the high temperatures and pressures of the mantle. For the reader unfamiliar with thermal and seismic equations of state, but who wants to delve further into this subject for its practical ramifications (e.g. Núñez-Valdez *et al.* 2013), understanding the interrelationship between elastic constants and thermodynamic variables is quite important.

### Thermal parameters

A key element of any geodynamical model is the thermally induced buoyancy. For our purposes, one element of this buoyancy, the thermal expansivity *α*_{th} ≡ −∂ ln*ρ*/∂*T* is a simple elastic parameter that can be derived by mineral physics (e.g. Anderson *et al.* 1992; Katsura *et al.* 2010) and it expresses the volume changes that occur due to heat content. The radial dependency of *α*_{th} can be adopted, for example, from Anderson *et al.* (1992) for the upper mantle and from Chopelas and Boehler (1992) for the lower mantle. Here we select the depth-dependence derived by Katsura *et al.* (2010). Our goal is to isolate anomalous heat content through its influence on shear-wave speed. Thus, the essential parameters are the thermal expansivity, *α*_{th}, and *c*_{1} or *γ*_{μ}. These parameters were the essence of the Ivins and Sammis (1995) method and implicitly enter in all of the approaches in this chapter.

To describe anomalous heat content it is important to have a 3D description of the temperature distribution in the mantle. We define anomalous temperature as *δT*(*r*, *θ*, *ϕ*), that part of the temperature that deviates from the globally based spherically symmetrical part. Stated formally: *δT*(*r*, *θ*, *ϕ*) ≡ *T* − *T*_{0}(*r*), wherein *T*_{0}(*r*) is the global spherically averaged temperature at radius *r*. This decomposition is standard procedure in computing free convection in planetary interiors (e.g. Durney 1968). To infer effective viscosity, all approaches require, implicitly or explicitly, a reference temperature profile in order to specify a radial variation of viscosity and a distinct laterally varying viscosity. Under the approximations developed in Ivins and Sammis (1995) these viscosity and temperature fields may be expanded in orthogonal functions. Hence, the relations for the viscosity and seismic velocity fields can both to be expanded in spherical harmonics and interrelated in the spectral domain. We do not, however, give the details of this procedure in this chapter. Table 2 gives the *T*_{0}(*r*) values used for all three approaches where we also show non-global averages at 150 km depth.

The elastic parameters can be combined to form the desired derivative that relates seismic velocity anomalies to temperature anomalies:

What follows is a discussion of the three approaches to deriving lateral viscosity variations from a seismic model. This requires an estimate of the global average temperature, *T*_{0}, at each depth. Three depths are chosen for the intercomparison: *z* = 150, 300 and 550 km. Because the reference seismic and viscosity models are global, it is essential that the reference temperatures represent global averages. The 150 km depth lies within the conducting lithosphere for some locations. Thus, for *T*_{0} at 150 km, we use the global average estimate of 1473 ± 40 K from Stacey and Davis (2008). Should non-global *T*_{0} values for relatively hot or cold subcontinental mantle be assumed (see Table 2), they will seriously bias the predictions of deviations from a purely radially dependent viscosity. We employ estimates of the isentropic, convecting, variation with depth for our value of *T*_{0} for the deeper depths, arriving at 1722 ± 42 K for 300 km depth and 1930 ± 50 K for 550 km (see figs 3 and 4 in Katsura *et al.* 2010). Karato and Karki (2001) pointed out the potential pitfalls of deriving seismic constraints on either mantle temperature and/or driving buoyancy in geodynamic models without *a priori* correction for anelasticity. We return to this issue in the section discussing Approach 2.

### Road map for the three approaches

Karato and Spetzler (1990) and Karato (1993) demonstrated that all realistic mantle rock mineral assemblages will have elastic-wave velocities diminished by thermally activated defect processes, so-called anelastic dispersion (e.g. Jackson *et al.* 2002). This must be accounted for in seismic studies, and the three approaches generally do so in different ways. Some elements of the approaches are summarized in Tables 1 and 3. Approach 1 uses a method that was partly described in the maps of the lateral variation in mantle viscosity derived using a combination of the scaling described by Ivins and Sammis (1995) with the modifications described by Wu *et al.* (2013). As a consequence, a number of features of Approach 1 and Approach 2 are similar, and we shall attempt not to be redundant in describing these overlapping features. Each method does apply an anelastic correction, as first suggested by Karato (1993). The latter is a necessary feature of deriving variable temperatures from seismic tomography. However, the detailed application of this correction differs between approaches. For Approach 1 and Approach 3, the anelastic correction simply follows from the table recommended by Karato (2008, table 20.2, p. 376). For Approach 2 we assemble the components of the anelastic correction from the general relationships for diffusion-controlled energy loss as described by Karato and Karki (2001) and in the monograph of Karato (2008, equations 18.5, 18.6, 20.12 and 20.13). We also rederived the asymptotic relationships of Karato (2008, equation 20.15b and discussion). This revealed that the asymptotics come with the added caveat that the frequency-dependent power-law exponent (*α*) must be restricted to *α* ≤ 0.3. Without the asymptotics, the root of the anelastic theory and its temperature dependence is essentially identical to that employed by Guerri *et al.* (2016).

For our intercomparison exercise we generally assume rather different flow laws to be operative. Approach 1 assumes diffusion creep as described and quantified by Hirth and Kohlstedt (2004), and Approach 2 assumes GBS dislocation creep as described and quantified by Hansen *et al.* (2011). Readers that may be interested in the derivation that leads to the linear relations between dln*ν*_{S} and *δT*, and the viscosity variability that we describe below,

Approach 3 is unique in that it uses flow laws in a composite form, combining olivine dislocation, diffusion and GBS micromechanics for calculations at 150 and 300 km depths, and similarly for ringwoodite at 550 km for dislocation and diffusion. Other differences that distinguish Approach 3 are summarized in Tables 1 and 3. Approach 3 uses thermoelastic plus anelastic scaling to *δ*ln*ν*_{S} that is identical to Approach 1.

### Inference from modelling geodetic data

The construction of trends from the time series of geodetic observations of the movement of the crust for the past three decades (e.g. Dietrich *et al.* 2004; Bevis *et al.* 2009) provides essential information for the Antarctic GIA modelling community. Inferences of mantle viscosity beneath Antarctica from geodetic observations are useful for fine-tuning, and even calibrating, the parameters used in mapping seismic velocity to viscosity. Unfortunately, these estimates generally require some assumption about ice-load history, which is not well understood, and thus different estimates for the same region can vary significantly. Nield *et al.* (2014) used GPS-measured uplift following the 2002 break-up of the Larsen B Ice Shelf to constrain the upper mantle viscosity beneath the northern Antarctic Peninsula to about 1 × 10^{18} Pa s or less. A more recent GIA analysis by Samrat *et al.* (2020) also favours a quite low shallow mantle viscosity in this region (0.3 × 10^{18}–3 × 10^{18} Pa s) and a modest upper mantle viscosity of 4 × 10^{20} Pa s deeper than 400 km depth. For the northern Antarctic Peninsula, Ivins *et al.* (2011) used a different GPS dataset, and Gravity Recovery and Climate Experiment (GRACE) data. The modelled GIA uplift also included an ice model with historical retreat (*c.* 100–150 years) and determined a best-fit mantle viscosity of 4 × 10^{19}–7 × 10^{19} Pa s. Bradley *et al.* (2015) and Wolstencroft *et al.* (2015) estimated upper mantle viscosity at 1 × 10^{20}–3 × 10^{20} Pa s across the southernmost Antarctic Peninsula and western Weddell Sea based on GPS measurements and improved ice-load history. Zhao *et al.* (2017) also modelled regional mantle viscosity in the central Peninsula and found values in excess of 2 × 10^{19} Pa s with an ice-loss history documented back to 1966 using remote-sensing data. Late Holocene beach uplift data was modelled by Simms *et al.* (2012) for the South Shetland Islands glacial advance and retreat, and associated viscoelastic GIA response. They determined a best-fit upper mantle viscosity value of 1.0 × 10^{18} Pa s, applicable to the asthenospheric mantle adjacent to the actively spreading Bransfield Strait Basin (e.g. Lawver *et al.* 1995). In summary, estimates of upper mantle viscosity beneath the Antarctic Peninsula range from less than 10^{18} to about 3 × 10^{20} Pa s, with lower values generally predominating in the northern and western areas.

A few viscosity estimates are also available for other areas of West Antarctica. Nield *et al.* (2016) suggested, based on observed responses to the stagnation of the Kamb Ice Stream, that upper mantle viscosity beneath the Central Siple Coast cannot be less than 1 × 10^{20} Pa s. Barletta *et al.* (2018) modelled extremely rapid GPS uplift rates in the Amundsen Sea Embayment region. After corrections for the elastic response from current ice-mass loss are applied, an upper mantle viscosity that is lower than the canonical upper mantle values (2–6 × 10^{20} Pa s) can be inferred. Their best-fitting model revealed a mantle viscosity of 4 × 10^{18} Pa s from 65 to 200 km depth, underlain by 1.6 × 10^{19} Pa s from 200to 400 km depth and 2.5 × 10^{19} Pa s from 400 to 660 km. Using a shorter ice-unloading history than Barletta *et al.* (2018), a GIA model proposed by Powell *et al.* (2020) that employed seismic scaling, suggests that viscosity as low as several times 10^{18} Pa s in the shallow upper mantle and several times 10^{19} Pa s down to a 900 km depth might be responsible for the rapid GPS uplift. Inevitably, these estimates of Antarctic mantle viscosity will improve as our constraints on past ice change become more robust (e.g. Nichols *et al.* 2019; Johnson *et al.* 2020).

We caution the reader that here we have discussed the ability of Maxwell Earth models, with a single parameter, *et al.* 2020; Adhikari *et al.* 2021).

In the following three sections we describe the methods and calibration for each of the three intercomparative approaches. Approach 2 also presents a disambiguation of the anelastic correction parameters. Here, analytical representation facilitates studying the sensitivity to all parameters. It is important to provide the reader with a concrete set of principles that the intercomparisons obey. Table 3 provides some additional rules that we follow. These supplement those presented in Table 1. The next three sections further detail those rules in the intercomparison that are defined in Table 3. Before we proceed, the background, or globally averaged effective viscosity at depth ^{20} Pa s, a value found for the upper mantle using the Antarctic GIA model of Ivins *et al.* (2013). Approach 2 tests the rheological law assumed at a depth of *z* = 225 km and compares it to the most reliable estimates that are available from the northern hemisphere. In Approach 2,

## Approach 1

### Methods for Approach 1

We now pursue computation of the ratio *M* ≡ log_{10}*e* and *a* ≡ (*E** + *PV**)/*R*_{G}. Here we recognize the exponent from the strain-rate equation (9) or the effective viscosity, *η*_{σ}, as in equation (14). Provided that one recognizes the qualifications for writing equation (14), the approximation (equation 25) is not particularly tied to diffusional creep laws, as it is equally valid for any dislocational creep law. It is assumed that there is an *r*, *θ*, *ϕ* dependence associated with log_{10}Δ*η* and *δν*_{S}, but only radial dependency with the other parameters. Clearly, the assumption of material homogeneity is implied in equations (9), (14) and (12), and is also implied in the practical application of the partial derivative relations derived from the elasticity relations for a realistic Earth model.

Wu *et al.* (2013) introduced a formulation that adds the anharmonic contribution of equation (24) to the contribution from anelastic dispersion:

This results in:
*β* discussed below. Karato (2008) provided approximations for the two derivatives on the right-hand side of equation (26) with values given in his table 20.2.

Computing viscosity using this method requires choices of global reference models for shear velocity, temperature and viscosity. Since the ANT-20 Antarctic seismic model was developed by Lloyd *et al.* (2020) from the global 3D model 362ANI, which was based on the 1D model STW105 (Kustowski *et al.* 2008), we use STW105 as the shear-velocity reference model. The temperatures are common to all three methods and are discussed in the previous section. Conductive temperature in the lithosphere is from a thermal boundary layer geotherm of Stacey and Davis (2008). The intersection of this profile and that of an adiabatic mantle occurs at about 200 km where we interpolate the geotherm onto the adiabatic temperature gradient of Katsura *et al.* (2010).

The computation requires rheological parameters for the construction of *a* in equation (25). For the activation energy, *E**, Hirth and Kohlstedt (2004) suggested a value of 375 kJ mol^{−1} for either wet or dry diffusion creep. While the activation volume, *V**, is poorly determined from experimental work, Hirth and Kohlstedt (2004) suggested that it could range from 2 × 10^{−6} to 10 × 10^{−6} m^{3} mol^{−1} for dry diffusion creep and *V** ≤ 20 × 10^{−6} m^{3} mol^{−1} for wet diffusion creep, so we use 5 × 10^{−6} m^{3} mol^{−1} in this calculation. These values are used for olivine mineralogy as it forms the dominant mineral phase above 410 km depth. At a depth of 550 km, we assume ringwoodite mineral assemblage, and hence laboratory experiments with this silicate phase are needed. Experiments by Fei *et al.* (2017) determined dislocation annealing-based mobility and its Arrhenius-dependence for both wet ringwoodite and bridgmanite at the appropriate mantle temperature and pressure. Here we use their estimates of activation enthalpy and the same value for activation volume as in Fei *et al.* (2017) to arrive at an estimate of the effective viscosity at 550 km depth.

The conversion of mantle shear velocity to mantle viscosity relies on the assumption that temperature variations are responsible for variations in both shear velocity and viscosity. However, other factors, particularly variations in mineralogy and volatile content, will also influence shear velocity. Here, the factor *β* ranges from 0 to 1, and is meant to account for effects on shear-wave velocity other than temperature (Wang *et al.* 2008). In the GIA analysis of Wu *et al.* (2013), *β* is used to relax the amplitude of the scale factor in fits to GIA-related datasets. In GIA analysis it forms a way of assessing the relative role of temperature in comparison to other factors affecting both the seismic velocity and log_{10}Δ*η*.

### Calibration of Approach 1

The recent Maxwell viscoelastic models for GIA developed during the past decade in Antarctica are used for calibration (see the ‘Inference from modelling geodetic data’ subsection earlier in this chapter). These model-based estimates are a potential source for placing realistic constraints on the parameter choices for the viscosity conversion. The extremely low upper mantle viscosity estimates of 1 × 10^{18}–4 × 10^{18} Pa s in the northern Antarctic Peninsula and the Amundsen Sea Embayment can only be approximately fitted using the IJ05-R2 upper mantle reference viscosity model of Ivins *et al.* (2013), which has the lowest continent-wide 1D upper mantle viscosity (2 × 10^{20} Pa s). Providing an acceptable match to these GIA-modelled observations also requires that shear-velocity variations result entirely from thermal variations (i.e. *β* = 1), at least in West Antarctica. A values of *β* close to 1 is compatible with evidence from compositional modelling which indicates that temperature variations have the greatest effect on mantle shear velocities, as realistic compositional variability produces only modest changes in seismic velocity (e.g. Goes *et al.* 2000; Lee 2003) and melt extraction produces only minor velocity increases in the upper mantle (Afonso and Schutt 2012).

A complete absence of compositional influence on shear velocities is unreasonable in East Antarctica, a craton that is underlain by thick continental lithosphere. On a global comparative basis, such lithosphere has been shown to be depleted, with high Mg/Fe ratios producing fast seismic velocity anomalies (e.g. Jordan 1981; van der Lee 2001; Lee 2003). Thus, we use *β* = 1, but the high mantle shear velocities in the upper reaches of the East Antarctic mantle are reduced by an appropriate factor prior to viscosity conversion. Assuming that depletion affects mantle shear velocity at the 2% level, Approach 1 makes a depletion correction by subtracting up to 1.7% from East Antarctic velocity anomalies greater than 3%, with the correction increasing with increasing velocity anomaly. The magnitude of this correction is based on both xenolith data from outside of Antarctica (Lee 2003) and petrological modelling (Afonso and Schutt 2012).

Clarification of the role of globally averaged *T*_{0} is important. Thermal convection in the upper mantle causes the globally averaged temperature gradient to deviate from an isentropic geotherm. However, using a common temperature for *T*_{0}(*r*) allows internally consistent results to be employed that come from high-pressure experimental and theoretical thermal elasticity. In Table 2 we give values for *T*_{0} that span relatively hot and cold subcontinental mantle environments and a global average estimate. This is important because in the shallow upper mantle, where the anelastic influence on *δ*ln*ν*_{S} is largest, the total log_{10}Δ*η* predicted is sensitive to the choice of geotherms, as shown in Figure 5 where we employ the ‘hot’ and ‘cold’ examples alongside the global average (see Table 2) as a function of the activation enthalpy for anelastic mobility.

## Approach 2

### Methods for Approach 2

A second method is now introduced that builds on the formulation of Ivins and Sammis (1995) as written in equation (25), but accounts for anelastic effects, as represented in the second term on the right-hand side of equation (26). Here we will disambiguate the components of the anharmonic (‘ah’) and anelastic (‘an’) contributions to the velocity anomaly, to aid in deciphering sources of error and bias.

The observed shear-wave velocity anomaly can be partitioned as:

All of the approaches essentially assume that

One might ask logically: to zeroth order, what is the temperature deviation from spherically averaged value that is caused by anharmonicity? A zeroth-order estimate of the thermal anomaly is:

We should note that derivations of the theoretical anharmonic estimates of *δT*, as in equation (29), depend on the exact flavour of the experimental and theoretical representations used, but at a basic level each must deal with the Grüneisen parameter and the temperature derivatives of the elastic moduli (e.g. Isaak *et al.* 1992; Núñez-Valdez *et al.* 2013).

The anelastic effect depends on the seismic shear-wave attenuation, *Q*_{S}. The physics describing this attenuation is more nuanced than is the description of seismic wave velocity using the theory of thermal elasticity. At the grain-boundary level and interior to rock grains, some strain during seismic wave propagation involves motions that cannot be described with elasticity because they involve microstructural imperfections. Although the causal mechanical interactions at grain boundaries and of interior dislocation distortions and point defect motions may be numerous, one feature is clear: they are highly sensitive to temperature through thermally activated processes (e.g. Schwarz and Granato 1975). Hence, these processes are also described by an Arrhenius temperature-dependence with an activation enthalpy,

Karato (1993) found an approximate solution for the amplitude of the incremental diminution in velocity that is a function of this temperature anomaly and which is related to the observed shear-wave seismic *Q*_{S} and the activation enthalpy:

Equation (30) shows the strongly non-linear relationship between seismic velocity and the thermal anomaly ratio *δT*/*T*_{0}. This particular form is strictly valid for the case of frequency-independent *Q*_{S}. Such a lack of dependency is consistent with the construction of ANT-20 using the Durek and Ekström (1996) mantle *Q*_{S} model. A rederivation of the frequency-dependent equivalent, however, shows that the approximation by equation (30) is nonetheless accurate to 10% provided that the power exponent of the frequency-dependence is less than 1/2, consistent with experimental estimates (Jackson and Faul 2010; Faul and Jackson 2015). We term this *et al.* (2014, equation 26, therein) for high-temperature background dissipation. This allows a first-order anelastic correction to be applied to the observed S-wave anomaly map. The anharmonic anomalies can be calculated from the observed S-wave anomalies (*δ*ln*ν*_{S}) as:

This result is essentially all that is needed to employ the method originally proposed by Ivins and Sammis (1995), which lacked corrections for temperature-dependent modulus dispersion (anelasticity). With this correction, we simply note that:
_{10}Δ*η*(*r*, *ϕ*, *θ*).

The advantage of computing the velocity anomaly explicitly with equation (31) is that the dependence on the seismic quality factor, *Q*_{S}, can be appropriately tailored to any attenuation model, whereas this cannot occur when employing the derivatives [∂ln*ν*_{S}/∂*T*]_{total} tabulated in Karato (2008). There is little fundamental difference in the way that Approach 1 and Approach 2 treat anelasticity. However, for Approach 1, all of the anelastic corrections are captured in the single term [∂ln*ν*_{S}/∂*T*]_{total} and these may be recovered from a table of values supplied by Karato (2008). For Approach 2, this term is disambiguated into its constituent parts, thus aiding in developing formal uncertainties as they relate to the formal uncertainties of those constituent parameters.

### Parameter sensitivity

A large body of experimental results quantify the spectrum of transient relaxation mechanics associated with thermally activated anelastic dispersion for olivine aggregates at high temperature (1173 ≤ *T* ≤ 1473 K) and pressure (0.2 GPa) (e.g. Jackson *et al.* 2002, 2004, 2014; Jackson and Faul 2010; Faul and Jackson 2015). Included are least-square fits to the experimental parameters on ensembles of samples strained under oscillatory conditions. These yield values for the activation energy for anelastic dispersion,

There is considerable sensitivity of the computed values of log_{10}Δ*η* at our depths of interest to *c*_{1} and *δν*_{S} seismic anomaly. The grey zone represents the breadth of values recommended by Karato (2008) and which satisfy laboratory values reported by Jackson and Faul (2010). The vertical dash-dot lines in Figures 6 show a value of *et al.* (2014) using an alternative parameterization of torsional experimental data and modelling of deformation in dry olivine. The latter quantified _{10}Δ*η* in terms of the observed lateral variability in S-wave velocity. A logical deduction from Figures 6 and 7 is that progress in theoretical and experimental interpretations of seismic anelasticity does impact our ability to make predictions of log_{10}Δ*η* in Antarctica.

### Calibration of Approach 2

The calibration of Approach 2 uses only past inferences of effective upper mantle viscosity in the northern hemisphere. The quality and abundance of the RSL and geodetic data in Fennoscandia makes it an ideal location for calibration, considering that our reference parameterization is global in nature. Studies that rely on the exponential decay curves of Holocene near-field uplift are especially important in Fennoscandia and Hudson Bay (e.g. Mitrovica 1996; Kuchar *et al.* 2019). The strongest sensitivity of these studies to the effective mantle viscosity is in the mid–upper mantle depth range. An ensemble of Fennoscandian ice-sheet reconstructions by Nordman *et al.* (2015) demonstrated that the GIA model solutions for upper mantle viscosity using the sea-level emergence record along the river valley of Ångermanälven on the western coast of the Gulf of Bothnia are decoupled from uncertainty in the ice-loading model. Such decoupling is critical to recovering tight bounds on the upper mantle viscosity, since the exponential decay of the rate of land emergence can be tightly coupled to the viscosity (e.g. Klemann and Wolf 2005). This fact motivates calibrating the *η*_{0}(*r*) to the experimentally developed flow law in Approach 2.

We place the activation enthalpies from experiments by Hansen *et al.* (2011) that determined the flow law for GBS dislocation mechanism into the approximation for effective viscosity using equation (14), resulting in:

In Hansen *et al.* (2011), experimental data were fitted with *n*_{GBS} = 2.9 ± 0.3, *p*_{GBS} = 0.7 ± 0.1 and

At an average background deviatoric flow stress *σ* = 0.5 MPa and a mantle depth of 220 km, where *T*_{0} ≃ 1500 K (Davies *et al.* 2012), we can obtain quite plausible values for *V** = 18 × 10^{−6}m^{3} mol^{−1} , suggested by Hansen *et al.* (2011) is also employed. At a depth of 220 km, hydrostatic pressure is roughly 7.1 GPa and accounting for the confining pressures of the experiments (300 MPa), the effective activation energy can be determined to be *E** ≃ 439.6 kJ mol^{−1}. We find with a grain size, *d*, of 4.5 mm, as suggested for application to the mantle by Hansen *et al.* (2011), and the appropriate hydrostatic pressure at this depth, then *et al.* 2006; Nordman *et al.* 2015), sites of relatively cold lithosphere. We note that GBS rheology is not assumed in Approach 1, so comparisons between the predictions of Approach 1 and Approach 2 may reveal biases that are tied specifically to the assumption of flow mechanism.

To derive an estimate of the effective viscosity at 550 km depth, we use the parameters derived by Fei *et al.* (2017) for hydrous ringwoodite. Using the activation enthalpy determined for ringwoodite (and activation volume, temperature and pressure estimate, see Table 4), we can estimate viscosity only by somewhat arbitrarily selecting additional parameters that have been derived in experiments for olivine, specifically using equation (12) with *A* = 1600, *d* = 25.5 mm, *p* = 3, *α* = *ϕ* = 0 and *n* = 3.5, as *σ* = 0.2 MPa (cf. Hirth and Kohlstedt 2004). With these, we calculate *et al.* (2018). The estimates that we have made calibrate *H** and temperature at the depth estimations we employ for log_{10}Δ*η*. Ultimately, the predictions should also be tested against Antarctic GIA inferences, as in Approach 1.

### Uncertainty quantification

We give predictions of the formal uncertainties in lateral variations of viscosity in Figure 8. The uncertainties are constructed from a combination of equations (25) and (29–32), yielding:
*H**, *α*_{th}, *T*_{0} and *c*_{1} are the important contributors to the uncertainty of the viscosity variations, and that they follow a Gaussian distribution, we can calculate the uncertainty of log_{10}Δ*η* as:
*σ*_{x} representing any one of the error estimates of the corresponding parameters in the sum. For example:

We give parameters assumed for the maps in Figure 8 in Table 4. These are guided by the parameter studies conducted as in Figures 6 and 7 (see also Fig. 9). Note that the estimates provide uncertainties for the viscosity conversion, but do not take into account any uncertainties in the shear velocities in the seismic model ANT-20. Furthermore, larger uncertainties than we can quantify here may be associated with unknown variability in composition, grain size or water. Deep uncertainty is also associated with the arbitrary selection of the effective Newtonian viscosity

## Approach 3

### Methods for Approach 3

The roots of Approach 3 have been applied previously to Antarctic rebound by van der Wal *et al.* (2015). The temperature, *T*_{0} + *δT*, is obtained following the description in the subsection ‘Thermal parameters’ and then inserted into a fully robust olivine flow law. The flow law allows diffusion, intragranular dislocation creep and GBS creep mechanics to operate simultaneously, each contributing to the macroscopic mantle strain (e.g. Ranalli 2001; Barnhoorn *et al.* 2011; Kohlstedt and Hansen 2015).

Dislocation creep is non-linearly related to stress, as has already been discussed. Viscosity, therefore, also carries this non-linear dependence, as is readily deduced from equations (9) and (16). van der Wal *et al.* (2013) formulated the composite rheology with the intention of solving time-dependent GIA loading–unloading problems in which effective viscosity evolves in time, as does the deviatoric stress and strain rate. However, here we desire a simpler approach, one that would facilitate intercomparison with the predictions of Approaches 1 and 2. We provide estimates for a constant stress state. Using the seismic velocity anomaly model ANT-20, the relative change in viscosity is computed as in equations (26) and (27). The spherically averaged global background temperatures used are identical to Approach 1 and Approach 2, as described in Table 2. The derivatives of the seismic velocity anomaly to temperature anomaly are provided as a function of depth in Karato (2008), thus accounting for both the anelastic and anharmonic contributions to S-wave velocity perturbations.

Unique to Approach 3 is that we add this microscopic contribution to the strain rate resulting from diffusion, dislocation and GBS mobility in the following way:
*et al.* (2011). Similar relations apply to the dislocation and diffusional component of the strain rate. To complete the evaluation of effective viscosity, the resultant

As in Approach 1 and Approach 2, we assume that the ringwoodite phase dominates at 550 km depth. Here, however, we select experimental results published by Kawazoe *et al.* (2016). In Figure 9 we show the lateral variability in viscosity that is predicted if we assume the parameters that are derived from experiments with ringwoodite at 18 GPa by Kawazoe *et al.* (2016). The experiments of Fei *et al.* (2017) used an annealing technique to obtain activation enthalpy values for dislocation mobility at temperatures and pressures relevant to the 550 km depth. However, annealing experiments do not apply a controlled shear stress and, hence, are not capable of determining a flow law. By contrast, the deformation experiments by Kawazoe *et al.* (2016) allowed a dislocation creep law to be derived, but with rather large uncertainties on the activation enthalpy. For a fit to the data in which the power exponent *n* was fixed to 3, *H*^{*} = 345 ± 90 kJ mol^{−1}; and when *n* was simultaneously solved for using the data, *H** = 279 ± 105 kJ mol^{−1}, with *n* = 2.4 ± 0.7. In Figure 9 we see that for a 2% S-wave anomaly, a discrepancy of no more than about 20% arises from the two differing fits to the data of Kawazoe *et al.* (2016). At 550 km depth, it is also probable that diffusion at grain boundaries contributes, and Shimojuku *et al.* (2009) provided diffusion constants for Si and O atomic species with activation enthalpies at 402 and 246 kJ mol^{−1}, respectively. We can use the former slowest diffusing species as rate controlling for diffusion creep in ringwoodite and form an estimate for the composite strain rate limited to

Unlike Approaches 1 and 2, wherein the full forms of the experimental flow laws are used for calibration purposes, the map views generated for Approach 3 for depths of 150 and 300 km attempt to use explicitly the experimentally determined prefactor and grain size (Figs 10c and 11c). For the ringwoodite mineralogy, it is more difficult to generate prefactors directly from the experiments and there are no experiments that collectively reveal the grain-size dependence. In Figure 12c, stress was set to 0.002 MPa and a self-consistent prefactor was derived.

### Calibration of Approach 3

The largest uncertainties in obtaining effective viscosity from Approach 3 comes from grain size, water content and stress (Hirth and Kohlstedt 2004; O'Donnell *et al.* 2017). Effective viscosity is very sensitive to these values because no average viscosity is assumed in the calculation of effective viscosity. To establish that the approach is realistically retrieving effective viscosity, model predictions can be compared to GIA data. The general composite approach with both linear (diffusion control) and non-linear (dislocation control) constitutive relations is employed in computing predictions of vertical land motion and compared to GPS observations in Antarctica (van der Wal *et al.* 2015), Fennoscandian RSL and GPS crustal motion data (e.g. Tushingham and Peltier 1992; Lidberg *et al.* 2010), and GRACE gravity trends (e.g. van der Wal *et al.* 2013). The computations are also used to estimate the GIA correction for the present-day mass balance of the Antarctic and Greenland ice sheets (e.g. Bouman *et al.* 2014; Xu *et al.* 2016). In general, a rheology with grain sizes of 4 mm or larger resulted in better fits to the uplift and gravity trend datasets. In contrast, a lower viscosity associated with smaller grain size and wet rheology predicts GIA signatures that underestimate the observed uplift rates. The dry rheology with 4 mm grain size is the preferred model in van der Wal *et al.* (2013) and is used here as the reference model. A stress of 0.5 MPa is assumed, which is a representative stress level for GIA-induced stress beneath the lithosphere locations underneath and near the margin of an ice sheet (van der Wal *et al.* 2010). Mantle stress levels, however, depend not only on the stress induced by GIA, but also on mantle convection (Bredow *et al.* 2021, this volume). There are some important caveats regarding our preferred values of grain size and water content. First, it is cautioned that these results are obtained with ice-loading histories that are implicitly based on 1D Earth profiles of linear viscosity. Progress with ice histories based on composite rheology is under investigation (Huang *et al.* 2019). Second, on the basis of xenolithic rocks recovered in Antarctica (see Chatzaras and Kruckenberg 2021, this volume), small grain sizes and wet olivine cannot be excluded (e.g. van der Wal *et al.* 2015; Chatzaras *et al.* 2016).

## Map views predicted by the three approaches

Viscosity variation in log_{10} units in plan view is predicted using the three approaches at the three intercomparative upper mantle depths. The predictions are given in Figures 10, 11 and 12.

### Synopsis of results for Approach 1

The resulting predictive model from Approach 1 is generally in approximate agreement with the ice-load geodetic viscosity estimates mentioned earlier (see Figures 10a, 11a and 12a, which show the three depth predictions). For example, in the Amundsen Sea Embayment, the modelled effective viscosity for the asthenosphere shallower than 200 km is 5 × 10^{18} Pa s, compared to 4 × 10^{18} Pa s estimated by Barletta *et al.* (2018). The modelled values for the rest of the upper mantle are somewhat higher than in Barletta *et al.* (2018), but still within a half of an order of magnitude. We note that in a GIA-LV model, Powell *et al.* (2020) inferred a viscosity in this same region of 2 × 10^{18} Pa s. The prediction of Approach 1 is unable to match the extremely low viscosity (*c.* 1 × 10^{18} Pa s) estimated for the northern Antarctic Peninsula and near Graham Land, by Simms *et al.* (2012) and Nield *et al.* (2014), yielding 7 × 10^{18} Pa s for 65–200 km depth,and higher viscosity values at greater upper mantle depths. At a depth of 550 km, lateral variability is further reduced, and predicted effective viscosity is close to 10^{20} Pa s. At this depth two regions, one in West Antarctica and the other in East Antarctica beneath the Wilkes Subglacial Basin, have a viscosity predicted to be close to 4 × 10^{19}–8 × 10^{19} Pa s.

### Synopsis of results for Approach 2

Figure 10b shows the predictions for Approach 2 at 150 km depth. Here much of East Antarctica has an effective viscosity, *η*_{σ}, of slightly above 10^{21} Pa s, rising to just above 10^{22} Pa s in a cratonic core region that extends 1100 km from the South Pole into the interior of East Antarctica and along the Wilkes Subglacial Basin. At 150 km beneath West Antarctica, *η*_{σ} is predicted to be between about 2 × 10^{19} and 10^{20} Pa s. A viscosity as low as 10^{18} Pa s is predicted SW of the Balleny Trough in the Indian Ocean and 20° latitude north of Marie Byrd Land in the Pacific Ocean sector, each being far from the Antarctic continent. The 1−*σ* uncertainties propagated from the experimentally based parameters might allow for the West Antarctic mantle effective viscosity to be lowered to just below 10^{19} Pa s.

At 300 km depth (Fig. 11b) the East Antarctic pattern is similar to 150 km. However, at this depth a large swath of the subcontinent extending from the coast adjacent to Wilkes Subglacial Basin to the southern Antarctic Peninsula and Weddell Sea region is essentially Fennoscandian in character (*η*_{σ}) takes on values near 3× 10^{20} –8 × 10^{20} Pa s. A portion of coastal West Antarctica and Marie Byrd Land may be in the range 0.9–10 × 10^{19} Pa s, not too different in value from 150 km depth prediction, but now of smaller overall regional extent. Again, the 1−*σ* uncertainties could allow predictions to descend to values slightly below 10^{19} Pa s in a small section of the mantle beneath Marie Byrd Land. (At this 300 km depth beneath Marie Byrd Land, P-wave tomography by Lucas *et al.* (2020) shows this region has a robust slow anomaly relative to surrounding West Antarctic mantle). The variability predicted at 550 km depth (Fig. 12b) is much reduced, with viscosity everywhere below 10^{21} Pa s, and is similar to that predicted by Approach 1.

We anticipate for all approaches that there will be large viscosity contrasts located at the upper half of the upper mantle, for it is here where convective temperature and partial melts may express themselves strongly (e.g. Yuen *et al.* 1993; Vacher *et al.* 1996). Indeed, we see that about four orders of magnitude characterize the prediction shown for each part of Figure 10 at 150 km depth when considering both continent and oceanic regions in map view. Excluding variability beneath oceanic crust in Figure 10, the variability reaches approximately three orders of magnitude. The largest variation occurs between coastal Marie Byrd Land, a region known for active volcanism, and the East Antarctic craton. At 300 km depth a magnitude variability of four orders is also reached if we include oceanic mantle. Hot oceanic upper mantle may be the site of highly channelized flow, at least in some interpretations of *δ*ln*ν*_{S} at this depth (e.g. French *et al.* 2013). Our focus, however, is on mantle beneath the continent of Antarctica where it is possible to retrieve corroborative viscosity information using bedrock geodesy.

At a depth of 550 km we predict that lateral variability is much reduced and, in fact, seems characterized by 1–1.5 orders of magnitude. This is caused by the smaller activation enthalpy for ringwoodite (see Table 4). We cautiously note that while Fei *et al.* (2017) discovered substantial variability, possibly from water content, we have not modelled this possible origin for lateral variability in rheology, for seismic models like ANT-20 are unlikely to deliver the required quantitative information.

Confining the assessment to West Antarctic mantle, the contrast in effective viscosity reduces to a little less than three orders at 150 km and 300 km depths. We carefully note that examination of Figure 6a reveals that a reduction in the anelastic activation enthalpy, ^{−1}, assumed for the predictive maps, to 307 kJ mol^{−1}, a value preferred in the summary of Jackson *et al.* (2014), increases the log_{10} variability by a factor of 1.56 at 150 km depth. Regardless of which value is assumed, Approach 2 predicts strong contrasts that will influence various ice history–geodetic solutions for mantle viscosity estimated through comparison to GPS crustal motions. This is because the regional viscosity has an exponential control on the isostatic decay time. Note that the preliminary 1−*σ* errors at 150, 300 and 550 km do not alter this conclusion.

### Synopsis of results for Approach 3

The predictions using Approach 3 are distinctive in that they are constructed with the assumption of composite micromechanics and explicit adoption of the corresponding flow laws, at least at depths of 150 and 300 km. At a depth of 150 km the cratonic mantle beneath East Antarctica is predicted to be 10^{22}–10^{23} Pa s, and possibly higher in some locations (Fig. 10c). This high viscosity extends well into the surrounding oceanic environment. A small portion of the Antarctic Peninsula and Marie Byrd Land is predicted to be near 10^{19} Pa s, and for a small region well offshore beneath the Balleny Trough the effective viscosity is predicted to be a few times 10^{18} Pa s. At this same depth there is a three orders of magnitude gradient in viscosity predicted on-continent between the southern Transantarctic Mountains and Marie Byrd Land.

Viscosity at 300 km depth is predicted to be a few times 10^{18} Pa s in oceanic regions (Fig. 11c). This prediction of sub-oceanic mantle viscosity is a ubiquitous feature, except in the Weddell Sea region. The Marie Byrd Land viscosity prediction is close to 5 × 10^{18} Pa s. A narrow corridor of viscosity beneath West Antarctica extends from northern Victoria Land to Thurston Island and west of the Bellingshausen Sea that is below about 6 × 10^{19} Pa s. A somewhat parallel band of Fennoscandian-like viscosity extends from just south of the Balleny Trough to the western Weddell Sea and northern Antarctic Peninsula. Also, East Antarctica is predicted to have a distinctively Fennoscandian character from the coast to 500–1000 km into the interior from 0° E to 90° E. At 550 km depth, Approach 3 predicts two pockets of relatively low viscosity (≥3 × 10^{19} Pa s) beneath Marie Byrd Land and Wilkes Subglacial Basin (Fig. 12c). Nearly all other regions are near Fennoscandian values, except for small pockets just west of the Syowa Coast and in the Antarctic Peninsula region where viscosity values are predicted to be in the range 10^{21}–10^{22} Pa s range.

## ANT-20-based viscosity v. model inferences

It is useful to display the GIA model results and the predictions of scaled effective viscosity at 150 and 550 km depths alongside one another, as is done in Table 6. Here we have not listed GIA models and our predictions at 300 km depth because the cases of 1D model viscosity at that depth are identical to the value assumed at 150 km. In Table 6, only the study by Powell *et al.* (2020) assumed a GIA-LV (3D) model. Differences between GIA model predictions can arise from differences in the details of the GPS time series employed, the ice-load history or depth parameterization. Given the underlying deep uncertainty in parameterizing the scaling factor *and* the fact that the values of the viscosity that best fit geodetic data depends on the selection of the model thickness of the uppermost mantle layer (e.g. Hu and Freymueller 2019), we caution against overinterpreting Table 6.

### Assessment using ice history–geodetic estimates of viscosity

Given the simplicity of our assumption that seismic tomography can be used to infer temperature, and therefore an effective viscosity, it should not be too surprising to find deficiencies in correlating mapped 1D viscosity to mapped shear-wave speed. Lateral heterogeneity in water or melt content, chemistry or anisotropy may also influence the viscosity and seismic velocity fields, and we have not attempted to account for these factors. However, when comparing our results from the three approaches we find general agreement with geodetic and palaeosea-level inferences of effective viscosity from GIA modelling. Basically, all three methods predict high viscosity beneath East Antarctic cratonic lithosphere and low viscosity beneath West Antarctica, with the centres of predicted low viscosity beneath Ross Island and Marie Byrd Land, sites of extensive Quaternary volcanism (e.g. Kyle 1990; Wilch *et al.* 1999; Smellie *et al.* 2021).

Low-viscosity regions are also imaged along the Amundsen Sea Coast and the Antarctic Peninsula. Furthermore, there are vast and continuous swaths of mid-level upper mantle where Fennoscandian-like viscosity (e.g. Mitrovica 1996; Lidberg *et al.* 2010; Nordman *et al.* 2015) is predicted, consistent with inferences of a number of continent-wide 1D GIA models (e.g. Whitehouse *et al.* 2012; Ivins *et al.* 2013; Argus *et al.* 2014). Among the most notable failures of the continent-wide 1D GIA models is their inability to predict large present-day vertical uplift on bedrock sites in the Amundsen Sea sector (e.g. Groh *et al.* 2012), several hundred kilometres from the Voigt-averaged S-wave anomaly minimum of ANT-20 beneath Marie Byrd Land. A low-viscosity region is predicted in this general region by each of the three approaches at 150 and 300 km. This is quite generally in accord with inferences of low viscosity that come from regional GIA modelling (e.g. Barletta *et al.* 2018; Kachuck *et al.* 2020; Powell *et al.* 2020). As we have discussed in reference to each of the predictions from the three approaches, there is discord on the magnitude of the inferred viscosity minimum. We do not speculate further on the causes of such a discrepancy. The list might include seismology, poor ice reconstruction, insufficient layering, low-temperature plasticity, transient creep or subtle effects of the heterogeneity itself. However, a positive result of our study is that the viscosity contrasts inferred from scaling seismic tomography to temperature, and then to viscosity do correlate well with the patterns of upper mantle viscosity inferred from GIA models that are well rooted in geodetic observation.

### Scrutiny of the Antarctic Peninsula

The northern Antarctic Peninsula (Graham Land) is one region that rather poorly correlates with the ANT-20 scaling to viscosity in at least two of the three approaches, as each are limited to treating lateral variability as dominated only by thermal control. The most optimistic correlation is predicted by Approach 1 at *z* = 150 km (see Fig. 10a). Models by Simms *et al.* (2012) and Nield *et al.* (2014) infer shallow upper mantle viscosity in the range of 0.7 × 10^{18} –2 × 10^{18} Pa s. The viscosity we present in this chapter based upon scaling to ANT-20 is nearly an order of magnitude higher than this at all depths we compute and for all approaches, presenting a challenge for interpretation.

In dealing with this enigma, three candidate causes are of primary consideration. First, there are potential changes in mantle grain size and volatile content caused by a youthful ridge subduction tectonics that characterize the northernmost Peninsula upper mantle environment (e.g. Bercovici *et al.* 2015). Secondly, is a resolution issue: a sharp variation in viscosity near the Bransfield Strait spreading centre might be captured in the geodetic data, but smeared out in the Voigt-average ANT-20 model. Generally, the model has approximately 150–200 km lateral resolution and about 50 km depth resolution at 150 km depth, but this resolution diminishes with depth. Thirdly is the possibility that both the GIA forward model set-ups and the experimental flow laws assume ‘steady-state’ properties, the former from the Maxwell viscoelastic assumption and the latter from creep experiments that have reached steady state, and therefore lack any transient behaviour. We are not yet in a position to distinguish between these three possibilities, and any further elaboration should be tempered by the fact that a higher-viscosity solution was derived from geodetic data by Ivins *et al.* (2011). The latter are seemingly compatible with two of the three predictions at 150 and 300 km depth, as these are between 1 × 10^{19} and 9 × 10^{19} Pa s (cf. Figs 10a, b & 11a–c). Additional seismic imaging of the upper mantle in this complex region might shed further light on these issues (e.g. Park *et al.* 2012).

## Discussion and conclusions

The three intercomparative approaches that we have used to convert the ANT-20 S-wave tomography of the Antarctic upper mantle to viscosity has been an exercise in testing for coherency in the resulting maps of log_{10}Δ*η*(*r*, *θ*, *ϕ*), noting the considerable differences in the underlying assumptions and procedures. The contrasts in the approaches are substantial, as can be deduced from the summary Tables 1 and 3. In fact, what we have sampled are methods that are in current use, and therefore give the reader a sense of the spread in results that underlie any transfer function for S-wave images to log_{10}Δ*η*(*r*, *θ*, *ϕ*) using laboratory-based steady-state flow laws. A more systematic study might be designed in the future that would consider a somewhat more broadened scope than that which we have developed here.

A key part of this chapter on rheology has been to revisit the issue of anelasticity in deriving quantitative relationships between lateral variation in effective mantle viscosity and S-wave velocity anomalies. Central to the physics of anelasticity are the grain-scale deformation processes associated with seismic wave motion. During the past two decades, torsional oscillation studies of GBS and static creep deformation in olivine have revealed that diffusionally assisted interionic mobility is a key feature of seismic attenuation (e.g. Jackson 2000; Marquardt and Faul 2018). Here we have taken the basic approach outlined in the monograph of Karato (2008). As the process of seismic wave attenuation involves grain-boundary assisted diffusion, an activation enthalpy of a transient creep, *et al.* 2002, 2014; Jackson and Faul 2010; Faul and Jackson 2015). However, these experiments are performed under an array of conditions (e.g. grain sizes, pressures, temperatures and synthetic preparation procedures) and, hence, there is a relatively broad range of

Our results reveal that among the competing parameters for calculations at *z* = 150 and 300 km for a 2% S-wave anomaly, _{10}Δ*η* that is predicted. Approaches 1 and 3 have selected values for the anelastic correction that are taken directly from the numerical values as a function of depth as supplied by Karato (2008). Karato (2008) cautioned about the sensitivity to the background *Q*_{S} model. However, we find that these cause negligible differences between the predictions from the three different approaches for Antarctic upper mantle. The background *Q*_{S} model used in construction of ANT-20 (*Q*_{S} = 70) differs only modestly from the model presented by Karato (2008) with *Q*_{S} = 80 at *z* = 150 km. In Figure 13 we explore the prediction for *δ*ln*ν*_{S} = −6%, a value characterizing the minimum S-wave velocity anomaly in ANT-20. Similar comparisons at depths of 300 and 550 km led to the same conclusions: that the differences in these particular *Q*_{S} models (𝒪 10%) play a negligible role. However, if the background *Q*_{S} model is larger by a factor of 2, then the inverted S-wave anomaly values increase (e.g. Lloyd *et al.* 2020, fig. S7). In addition, the influence on the anelastic correction would then be smaller, resulting in a larger prediction of log_{10}Δ*η*(*r*, *θ*, *ϕ*) associated with a thermal perturbation *δT*(*r*, *θ*, *ϕ*).

In Figure 13 we relax the rule of thumb suggested by Karato (2008) (0.8 ≤ *ξ* ≤ 1) in the relation *et al.* (2014) suggested a relaxation in the lower bound on *ξ* to allow *ξ* ≈ 0.5. Such lowering of *ξ* is also consistent with the analysis of McCarthy *et al.* (2011). A reduction to this value does have some consequence. For example, the map of Figure 10b shows viscosity variability weaker than that allowed in Figure 13. At _{10}Δ*η* would be raised by about a half an order of magnitude with respect to values in Figure 10b if using this lower activation enthalpy for thermally activated anelasticity. Potentially, ambiguities in the range of acceptable *a priori* lateral viscosity models derived from seismology.

One method that can account for the different creep mechanisms that correspond to laboratory-derived flow laws is to abandon the use of viscosity as a single controlling parameter, and employ the creep law of equation (9) directly into the initial-value boundary-value problem (IVBVP) code structure that solves geodynamic and GIA model set-ups. This, in fact, has been done in GIA models (e.g. Wu 1999; Zhong *et al.* 2003). Diffusion, dislocation and GBS creep can operate simultaneously in the mantle (Ranalli 2001; Kohlstedt and Hansen 2015), and a composite rheology has been used in mantle convection (e.g. van den Berg *et al.* 1995; Becker 2006; Dannberg *et al.* 2017) and in GIA models (e.g. van der Wal *et al.* 2013; Huang *et al.* 2019). In this case, a mean effective viscosity could only be defined as a retrospective summary of the class of full numerical simulations, should the fully non-linear constitutive equations be employed. As our goal has been to discuss the rheology of the mantle beneath the Antarctic ice sheet that can be derived from a comprehensive S-wave tomographic image, we must acknowledge this limitation. Tractable relations between S-wave anomaly and the single Newtonian-like effective viscosity make bold approximations in order to be easily compared to viscosity derived from ice history–geodetic reconstructions of GIA (e.g. Ivins *et al.* 2011; Nield *et al.* 2014; Barletta *et al.* 2018; Powell *et al.* 2020). The latter models are nearly all based on a Maxwell constitutive assumption and report a Newtonian viscosity value. Hence, the scaling methods explored in this paper have a practical utilitarian value since the preponderance of GIA models are 1D and assume a Newtonian viscous element in a Maxwell rheology. An inherent goal of developing coherent transfer functions between tomography and viscosity is to develop more advanced 3D structures for use in GIA models. For example, in a study of the Amundsen Sea Embayment, Powell *et al.* (2020) showed that the sensitivity of horizontal GPS station motions to 3D structure is substantial. For regions peripheral to the locus of ice-mass change, 3D structure can also influence vertical motions since bulge migration involves lateral transport of mantle material (Klemann *et al.* 2007).

We conclude that the basic philosophy of deriving scaling relationships is a coherent and rationale one. With quantification and propagation of laboratory-derived and seismically derived uncertainty, we are also able to describe some of the uncertainties in the scaling to effective regional viscosity. A future challenge may be to explore more fully various poorly understood biases, and to understand better the grain-size dependency and impurity effects in the deeper upper mantle for wadsleyite and ringwoodite mineralogy. The issue of lower mantle viscosity just below the 660 km discontinuity has important implications for the space gravimetric interpretation of present-day ice-mass changes in East Antarctica, as shown recently by Caron and Ivins (2020). The coming decade will see the development of increasingly sophisticated numerical models for GIA-LV, possibly even those that will employ adjoint methods (e.g. Crawford *et al.* 2018). Critical to these advanced 3D approaches is having high confidence in a starting model that is well rooted spatially. For this purpose, seismic mapping in the mantle will be a valuable tool. Among a number of caveats we have not touched on, the most important is the *other* implication of the mobility processes associated with the activation parameter *et al.* 2020). Transient rheology could help to explain the high amplitudes of geodetically determined responses to decade timescale ice-unloading events in Antarctica. How to deal with this short timescale rheology will be a major future challenge.

## Acknowledgements

We dedicate this chapter to the memory of Orson L. Anderson and Donald L. Anderson. We thank the editor, Adam Martin, Pippa Whitehouse and two anonymous reviewers for quite thorough and helpful suggestions.

## Author contributions

**ERI**: conceptualization (equal), data curation (supporting), formal analysis (lead), funding acquisition (lead), investigation (lead), methodology (equal), project administration (lead), resources (equal), software (equal), supervision (lead), validation (lead), visualization (equal), writing – original draft (lead), writing – review & editing (equal); **WVDW**: conceptualization (equal), data curation (equal), formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), project administration (equal), resources (equal), software (equal), supervision (equal), validation (equal), visualization (supporting), writing – original draft (supporting), writing – review & editing (supporting); **DAW**: conceptualization (equal), data curation (equal), formal analysis (equal), funding acquisition (equal), investigation (equal), methodology (equal), project administration (equal), resources (equal), software (equal), supervision (equal), validation (equal), visualization (equal), writing – original draft (equal), writing – review & editing (equal); **AJL**: conceptualization (supporting), data curation (equal), formal analysis (supporting), funding acquisition (supporting), investigation (supporting), methodology (supporting), project administration (supporting), resources (equal), software (supporting), supervision (supporting), validation (supporting), visualization (supporting), writing – original draft (supporting), writing – review & editing (supporting); **LC**: conceptualization (supporting), data curation (supporting), formal analysis (supporting), funding acquisition (supporting), investigation (supporting), methodology (supporting), project administration (supporting), resources (supporting), software (supporting), supervision (supporting), validation (supporting), visualization (supporting), writing – original draft (supporting), writing – review & editing (supporting)

## Funding

This research was carried out at the Jet Propulsion Laboratory (JPL), California Institute of Technology, under a contract with the National Aeronautics and Space Administration (NASA) and funded through NASA Post-doctoral Program, and through NASA's Earth Surface and Interior Focus Area Program (grant NNH15ZDA001N-ESI to E.R. Ivins), the Sea-Level Change Science Team (grant NNH16ZDA001N-SLCT to E.R. Ivins), the GRACE Follow-On Science Team (grant NNH19ZDA001N-GRACEFO to E.R. Ivins) and the NASA Cryosphere Program (grant NNH13ZDA001N-SLR to E.R. Ivins).

## Data availability

The data that support the findings of this study are available from Jet Propulsion Laboratory, California Institute of Technology, but restrictions apply to the availability of these data, which were used under licence for the current study, and so are not publicly available. Data are, however, available from the authors upon reasonable request and with permission of Jet Propulsion Laboratory, California Institute of Technology.

## Appendix A: A linearization of the thermal perturbation

In this Appendix we show the level of bias that is inherent to a linearization of *δT* in equation (7). The approximation is technically unnecessary, but it eases a direct use of spherical harmonics for simultaneously treating global tomography models and semi-analytic formulations of the GIA-LV problem (e.g. Martinec 2000; Tromp and Mitrovica 2000; Tanaka *et al.* 2011). Spherical harmonic representation of the field log_{10}Δ*η*(*r*, *θ*, *ϕ*) is then smoothed with the same resolution and fidelity as the seismic tomographic map of *et al.* 2005). Additional terms in equation (8) include

Comparison of the exact relationship (equation 7) and two levels of approximation are shown in Figure A1. Also shown here (dotted curve overlying the yellow dash-dots) is the method of Latychev *et al.* (2005), which is also a first-order approximation. The study by An *et al.* (2015) used a Rayleigh-wave dispersion analysis to determine shear-wave structure. It is of note that if we adjust the parameter *ɛ* of Latychev *et al.* (2005) in equation (4) to a value of 0.067 K^{−1} the predictions become identical to the first-order approximation that employs physics into the parameterization. This contrasts to the ɛ value assumed in the global modelling by Austermann *et al.* (2013), *ɛ* = 0.04 K^{−1}, and to one of the values used for Antarctica by Hay *et al.* (2017) (

In Figure A2 we show the positions from the steady-state thermal model computed by An *et al.* (2015) that were used in constructing Figure A1. Our method employed in this chapter uses a different set of relationships to compute temperature as a vehicle to estimate lateral viscosity contrasts. Nonetheless, we use Meijian An's model as an independent method for estimating *δT* and *T*_{0} for the Antarctic mantle. It therefore provides an independent way of testing the linear approximation for log_{10}Δ*η* widely assumed in the GIA-LV modelling community (e.g. Hay *et al.* 2017; Powell *et al.* 2020).

- © 2021 Jet Propulsion Laboratory, California Institute of Technology

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)