Sea levels and uplift rate from composite rheology in glacial isostatic adjustment modeling

Journal of Geodynamics 50 (2010) 38–48

Contents lists available at ScienceDirect

Journal of Geodynamics journal homepage: http://www.elsevier.com/locate/jog

Sea levels and uplift rate from composite rheology in glacial isostatic adjustment modeling Wouter van der Wal a,b,∗ , Patrick Wu c , Hansheng Wang d , Michael G. Sideris a a

Dept. of Geomatics Engineering, University of Calgary, 2500 University Drive, N.W., Alberta, Canada T2N1N4 DEOS, Faculty of Aerospace Engineering, Delft University of Technology, Delft, The Netherlands Dept. of Geoscience, University of Calgary, 2500 University Drive, N.W., Alberta, Canada T2N1N4 d Key Laboratory of Dynamic Geodesy, Institute of Geodesy & Geophysics, Chinese Academy of Sciences, Wuhan 430077, China b c

a r t i c l e

i n f o

Article history: Received 21 November 2008 Received in revised form 18 November 2009 Accepted 6 January 2010

Keywords: Glacial rebound Finite-element method Composite rheology Power-law creep Relative sea level Uplift rate

a b s t r a c t Since laboratory experiments point to the existence of both diffusion creep and power-law creep at realistic mantle conditions, a composite rheology in which both mechanisms operate at the same time so that their strains are combined might be more realistic and has in the past been shown to provide a better fit to observations of glacial isostatic adjustment (GIA) than purely linear rheology (Gasperini et al., 2004; Dal Forno et al., 2005). To further investigate the effect of such rheology on sea level curves and uplift rate resulting from GIA, composite rheology is implemented in the coupled Laplace-finite element method for a 3D spherical self-gravitating Earth. We vary the pre-stress exponent (assumed to be derived from a uni-axial stress experiment) and the Newtonian viscosity for a homogeneous mantle. Composite rheology is found to have a statistically significant better fit with observed relative sea level data than linear rheology (diffusion creep only) and non-linear rheology (dislocation creep only). For the best-fitting composite rheology model it is shown that in the mantle below the former ice sheet margin, stress is high enough for power-law creep to become dominant during melting and shortly thereafter, causing the model to behave mostly in a non-linear way. It is found that composite rheology, with the parameters investigated in this paper, not only provides a better fit to sea level data than nonlinear rheology but also slightly increases present-day uplift rate compared to non-linear rheology. This encourages application of composite rheology in GIA models that aim to improve knowledge of mantle rheology. Low uplift rates for composite rheology can be further increased by a large increase in ice thickness in North America at the expense of violating total melt-water constraints. A 1 or 2 ka delay in Laurentide glaciation and deglaciation increases uplift rates for all values of the pre-stress exponent investigated, while fit to a number of relative sea level observations in the Laurentide ice sheet is improved. Large increase in ice thickness disagrees with other observations (total melt constraints), therefore a delay in glaciation is a promising direction if global ice models are to be adjusted for a composite rheology. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction Most glacial isostatic adjustment (GIA) models are based on a linear relation between stress and strain rate, i.e., linear or Newtonian rheology, which has been successful in explaining a large number of GIA observations simultaneously, including relative sea level (RSL) data, crustal velocities, gravity rates of changes, polar wander and non-tidal acceleration of the Earth (e.g., Peltier, 1998, 2004). However, this success alone does not justify the use of linear rheology. Laboratory measurements show that dislocation creep,

∗ Corresponding author at: DEOS, Faculty of Aerospace Engineering, Kluyverweg 1, 2629 HS Delft, The Netherlands, Tel.: +31 15 278 2065. E-mail address: [email protected] (W. van der Wal). 0264-3707/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.jog.2010.01.006

which gives rise to non-linear or power-law creep, also operates in the mantle. Formally, power-law creep can be written in tensor form as (Ranalli, 1995): n−1

ε¯ ij = A  E

ij

(1)

where ε¯ ij is the deviatoric strain rate, ij is the deviatoric stress, =   n−1 E



(1/2)ij ij is the effective shear stress, n is the stress

exponent, A is a parameter determined from shear experiments and is a function of pressure, temperature and material properties. It is also useful to define the effective viscosity eff as (Wu, 1998): eff =

1 2A  n−1 E

=

1 3A∗ qn−1

(2)

W. van der Wal et al. / Journal of Geodynamics 50 (2010) 38–48

√ where q = 3E is the von Mises equivalent stress, and A*, the parameter determined from uni-axial experiments, is related to A by A∗ = (2/3(n+1)/2 )A (Ranalli, 1995). The parameter A* and the von Mises equivalent stress are introduced because the ABAQUS finite element model uses them as inputs for Earth deformation calculations, and previous works (Wu, 1995, 2001, 2002; Wu and Wang, 2008) have used A* instead of A. Laboratory experiments found that power-law creep (n > 1) operates when the stress level is high or the grain size is large (e.g. Goetze and Kohlstedt, 1973; Ranalli, 1991; Karato and Wu, 1993; Li et al., 2003; Cordier et al., 2004; Mainprice et al., 2005). In contrast, if the stress level is low or the grain size small, deformation occurs mainly through diffusion creep and the flow law is linear (n = 1). However, due to the large uncertainty of the conditions for the linear–non-linear transition (e.g. Ranalli, 2001), the presence of water in the mantle and post-perovskite in the lower mantle, microphysics alone cannot say definitively which part of the mantle behaves linearly or non-linearly. GIA observations have been used to address this question, see a recent review by Wu and Wang (2008) which we will summarize below. For a power-law mantle, Schmeling (1987) argued that GIA should see a linear creep law while mantle convection should see a non-linear creep flow. While this is true for RSL data near the center of rebound, Wu (1995) showed that RSL curves near the ice margin are diagnostic of non-linear rheology. This was confirmed in a 3D flat Earth model with a realistic ice history and ocean loading (Wu, 2001). The best fit to RSL data was obtained for an Earth model with power-law rheology restricted to the lower mantle (Wu, 2002; Wu and Wang, 2008). These studies limit non-linear rheology to certain layers in the mantle, while linear rheology is assumed in the other layers. However, there is no reason why linear and non-linear rheology cannot co-exist, and a so-called composite rheology seems a better approximation to real Earth deformation (Ranalli, 2001; Korenaga and Karato, 2008). A composite rheology can be constructed by the Ellis model (e.g., Bird et al., 1960), which sums the creep rate from linear and non-linear flow laws. It is used in mantle convection by, e.g., Parmentier et al. (1976) and a similar model is used by van den Berg et al. (1993) and Gasperini et al. (1992). The latter is the first application of composite rheology in GIA modeling and is followed by investigations in a series of recent papers. Gasperini et al. (2004) used a flat axisymmetric model and found that composite rheology was able to fit RSL data better than linear rheology. This is confirmed by Dal Forno et al. (2005) and Dal Forno and Gasperini (2007). Previous studies with composite rheology use a flat Earth geometry (Gasperini et al., 2004; Dal Forno et al., 2005; Dal Forno and Gasperini, 2007) while sphericity can become important for the Laurentide ice sheet especially for RSL data far from its center. Finally, self-gravitation and the self-consistent sea-level equation are neglected in previous studies (Giunchi and Spada, 2000; Gasperini et al., 2004; Dal Forno et al., 2005; Dal Forno and Gasperini, 2007). The sea-level can be important at the edge of the former ice sheet, where the ice attracted large amounts of sea water and where also linear and non-linear rheologies behave differently (Wu, 1995). Thus, it might be important to model self-gravitation when comparing different rheologies. In this paper, the effects of Earth’s sphericity, the self-consistent sea-level equation and selfgravitation are included in our study of composite rheology. A comprehensive comparison of linear, non-linear and composite rheology has not been shown in the previous studies. Dal Forno et al. (2005) and Dal Forno and Gasperini (2007) showed that composite rheology fits the RSL data significantly better than linear and non-linear rheology. However, a global misfit number does not illuminate the temporal and geographic spread of differences in sea

39

level curves, while the RSL behavior at a specific location is sometimes useful to distinguish between linear and non-linear rheology (e.g. Wu, 1995). Currently, it is not clear how composite rheology behaves for different values of A*, whether the predictions are closer to linear or non-linear rheology, and how predictions using a composite rheology depend on location and time. The answers to these questions can help us relate previous studies of non-linear rheology by Wu (2002) and Wu and Wang (2008) to studies of composite rheology (Gasperini et al., 2004; Giunchi and Spada, 2000; Dal Forno et al., 2005; Dal Forno and Gasperini, 2007). Moreover, present-day uplift rates from composite rheology have not been studied, while a known problem with non-linear rheology is the low uplift rates it provides (e.g. Wu, 1999). Because Peltier’s global ice models ICE-4G and ICE-5G are based on linear rheology they might require modification when using them in combination with a non-linear rheology which generally leads to faster relaxation during early deglacial time but lower present-day uplift rates. Wu (1998, 1999) showed that an increase in ice height improves fit with relative sea level data within the ice sheet margin. Wu and Wang (2008) found that a 2 ka delay deglaciation applied to ICE-4G makes for a better fit of sea level curves in the center of the former Laurentide ice sheet. However, it is not clear how modifications in the ice model affect the predictions of composite rheology. The answer to this question can contribute to an assessment of the performance of composite rheology and can also benefit future work aimed at making an ice model compatible with non-linear or composite rheology. In summary, this paper aims to address the following questions: (1) In a composite rheology, which regions of the mantle are dominated by non-linear flow and how does the dominance of non-linear rheology in these regions change in time as deglaciation proceeds? (2) How do sea level curves predicted by composite rheology differ from those predicted by non-linear and linear rheology? (3) Can a good fit to sea level data as well as present-day uplift rate data be obtained with composite rheology? (4) How can ice models be modified to improve fit of composite rheology with sea level data and present-day uplift rate in North America? In order to investigate these questions, we use the definition of composite rheology in which strain rates from diffusion and dislocation creep are summed, assuming that both processes can occur simultaneously: ε¯ ij =

ij 2

n−1

+ A  E

ij =

ij 2

+ 32 A∗ qn−1 ij

(3)

where  is the Newtonian viscosity, and other quantities are defined as in Eq. (1). This formulation is identical to Gasperini et al. (1992) and Giunchi and Spada (2000) except that we keep the creep parameter A (or rather A*) as input parameter instead of the transition stress as in these studies (see Dal Forno et al. (2005) for the relation between transition stress and the linear viscosity and the creep parameter A). We have also implemented a composite rheology in which the strain rates from diffusion and dislocation creep are not summed, but the strain rate of the medium is the larger strain rate of the two. However, we found that the differences between predictions from these two formulations are small, so we just use the formulation of Eq. (3) in this paper. For this study we make the assumption that the value of A (or A*) is the same throughout the mantle although in reality A (or A*) is dependent on activation energy, activation volume, temperature, water content, grain size and shear modulus (e.g., Turcotte and Schubert, 2002; Karato and Wu, 1993). Thus, our study should

40

W. van der Wal et al. / Journal of Geodynamics 50 (2010) 38–48

Table 1 Elastic parameters for the Earth modela . Layer

r (km)

Lith UM TZ LM1 LM2 Core

6371 6221 5971 5701 5200 3480

 (kg/m3 ) 3192 3442 3882 4527 5074 10,987

g0 (m/s2 )

E (×1011 Pa)



9.815 9.866 9.969 10.014 9.947 10.683

1.55 1.89 2.83 4.60 6.24 0

0.50 0.50 0.50 0.50 0.50 0

a r, the outer radius of the layer; , density; g0 , gravity at the top of the layer; E, Young’s modulus; , Poisson’s ratio.

be considered a preliminary study aimed to develop GIA models that are more in agreement with laboratory measurements and microphysics. Furthermore, we omit the background stress in this study because there is no clear knowledge of the location and the pattern of convection cells within the mantle. Moreover, Karato (1998) pointed out that the strain rate due to GIA is orders of magnitude smaller than that due to tectonics; consequently, the distance that dislocations move during GIA does not significantly exceed the average distance of dislocations. Thus, dislocation density is unlikely to change appreciably during GIA and there is little interaction between GIA induced stress and ambient tectonic stress. Gasperini et al. (1992), Giunchi and Spada (2000), Dal Forno et al. (2005) and Dal Forno and Gasperini (2007) include background stress as scalars. However, background stress is a tensor and can both increase and decrease the effective viscosity of power-law creep and hence the scalar sum of the background stress and the rebound stress only simulates the special case where the directions of the stress are parallel and of the same order of magnitude (Schmeling, 1987). We assume here that GIA sees steady state creep, following the argument of Ranalli (2001). Inversions of data on both long timescales (formation of the long-wavelength geoid) and the relatively short timescales of postglacial rebound result in the finding of a significant increase in viscosity from upper to lower mantle (although the viscosity contrast from GIA inversions is still under debate, see e.g. the recent review in Wolf et al., 2006). The agreement of inversions from GIA data with inversions of the longwavelength geoid argues against the importance of transient creep in postglacial rebound. 2. Finite element modeling with composite rheologies Our model couples the Laplace equation to a finite element method to solve for deformation on a spherically stratified selfgravitating incompressible Earth and oceans (Wu, 2004). In order to visualize easily the stress distribution during the glacial melting, we employ the simple axisymmetric model of Wu and van der Wal (2003) in Section 3. For the prediction of relative sea level at specific sites in Section 4 and uplift rates in Section 5, the 3D self-gravitating spherical Earth used by Wang and Wu (2006) is modified to include composite rheology. This 3D model is more realistic than previous studies of composite rheology because of the inclusion of both sphericity and self-gravitation, but it has the disadvantage of long computation time which reduces the range of models that can be investigated. To reduce computation time, the model has a spatial resolution of 2 × 2 degrees on the Earth’s surface. However, the coastlines and the ice margins are more winded than can be described by this coarse grid. Thus in the computation of relative sea levels, the location of one site (Newfoundland) had to be adjusted slightly (<1.5◦ ) to give it the proper relative position to the ocean and ice loads. A six-layer elastic Earth model is used to reduce computation time and memory required while still incorporating the major seismic discontinuities. Its elastic parameters are given in Table 1. It is the same as the model

used in Wu and Wang (2008) except that a 150 km thick lithosphere is used here instead of 115 km (in order to have a better agreement with the thick lithosphere in the North American craton, Peltier, 1984), and the model here is incompressible (Poisson’s ratio of 0.5). In this study, the creep parameters are taken to be constant in the mantle. For the linear rheology, the parameter A* in the mantle is taken to be 3.3 × 10−22 , 1.1 × 10−22 , 3.7 × 10−23 and 1.2 × 10−23 Pa−1 s−1 which correspond to a Newtonian viscosity of 1, 3, 9 and 27 × 1021 Pa s, respectively. For the non-linear rheology, A* (constant in the mantle) is varied between 3.3 × 10−33 , 3.3 × 10−34 , 3.3 × 10−35 and 3.3 × 10−36 Pa−3 s−1 , as this is found to be a reasonable range for power-law creep to match sea level data (Wu, 2002; Gasperini et al., 2004; Dal Forno and Gasperini, 2007; Wu and Wang, 2008). The stress exponent is taken to be n = 3 and is not varied, because Wu (2002) and Wu and Wang (2008) found that the sea level data cannot discern between n = 3 and n = 4. Thus, this paper considers 16 different composite rheology models, which are a combination of any one of the 4 linear rheologies with any one of the 4 values of A* for the non-linear rheology. For the axisymmetric model, the surface load consists of a 15◦ rectangular ice cap and complementary self-gravitating ocean. Ice thickness increases linearly during the 90 ka glaciation phase, then decreases linearly during the 9 ka deglaciation phase and is zero during the last 11 ka. For the 3D model, both the ICE-4G (Peltier, 1994) and ICE-5G (Peltier, 2004) ice loading histories are used to compute surface loads during deglaciation, but for the glaciation phase the load is assumed to increase linearly from the start of glaciation to the last glacial maximum (LGM). Note that in ICE-5G maximum ice heights are at 26 ka for some locations, with little change until 20 ka BP which is usually taken to be the LGM. The ice history in the axisymmetric model assumes the LGM at 20 ka before present (BP), more in accord with the ICE-4G model (Peltier, 1994). The self-consistent sea level equation is solved, accounting for ocean inflow into previously glaciated areas (Hudson Bay and the Gulf of Bothnia) as deglaciation proceeds. The ocean function is constant otherwise. Relative sea level observations from 128 sites from the Tushingham and Peltier (1991) database are used in misfit computations.

3. Spatiotemporal distribution of non-linear deformation The axisymmetric model is used to visualize which regions of the Earth model with composite rheology are dominated by linear and which by non-linear creep. Fig. 1 shows the elements below the elastic lithosphere in the axisymmetric composite model with A* = 3.3 × 10−35 Pa−3 s−1 and  = 3 × 1021 Pa s that behave as linear or non-linear creep. This model is selected as it is the best-fitting composite rheology model with respect to the relative sea level data (see Fig. 4 and discussion later). Three time steps are plotted: LGM (17 ka BP), end of melting (11 ka BP) and 3 ka after the end of the melting phase. The area that is dominated by non-linear creep grows from 17 ka BP, quickly expands in all mantle layers (but not below the center of the ice sheet) as deglaciation sets in, grows until the end of deglaciation and decreases quickly until all the mantle is dominated by linear creep after 7 ka BP. To understand this, we plot in Fig. 2 the time evolution of the von Mises stress at four locations ( = 0.25◦ , 10.25◦ , 15.25◦ and 30.25◦ from the center of the load) just below the lithosphere in Fig. 1. Note  from Eq. (3) that non-linear creep becomes dominant if qt = (1/ 3A∗ ) < q, where qt is the transitional von Mises stress. Fig. 3 plots qt versus log10 (A*) for four typical values of  in the mantle when n = 3 and for a range of values of A*. When the von Mises stress in the mantle is above these curves, non-linear creep dominates, otherwise linear creep dominates (see also Gasperini et al., 2004). In general, the peak von Mises stress level in the mantle

W. van der Wal et al. / Journal of Geodynamics 50 (2010) 38–48

41

Fig. 1. Area below the mantle layers of the axisymmetric model where linear rheology (light grey) or non-linear rheology (dark grey) dominates for a composite rheology with a Newtonian viscosity of 3 × 1021 Pa s and A* = 3.3 × 10−35 Pa−3 s−1 . The epochs shown are 17 ka BP or LGM, 11 ka BP or end of melting, and 8 ka BP. The rectangular ice cap extends from the North Pole to the solid radial line at 15◦ .

is typically of the order of a few MPa, so if the transitional von Mises stress is much higher than this value linear rheology will generally dominate. On the other hand, if the transitional von Mises stress is much lower than this value, non-linear creep will generally dominate. Fig. 3 shows that higher linear viscosity leads to lower qt and thus the deformation is more likely to be dominated by non-linear rheology. This is an important point that will help us understand the general behavior of more realistic models in the following sections. For A* = 3.3 × 10−35 Pa−3 s−1 and  = 3 × 1021 Pa s (Model I) it follows that qt ≈ 1.8 MPa; this value is indicated in Fig. 2 as the thick horizontal line. Therefore, when q evolves in time (Fig. 2) and rises above qt non-linear creep dominates; this is indicated as the dark grey region in Fig. 1. For example, in the load center just below the lithosphere the value of q is always below qt , thus the dominant creep is always linear (note that a large part of the 37 MPa total ice load at LGM is supported by the elastic lithosphere and that the von Mises stress is computed from deviatoric stresses). Fig. 1 also shows that the upper mantle below the ice center is dominated by linear rheology but the lower mantle is dominated by non-linear creep Fig. 3. Transitional von Mises stress qt versus log 10(A*) for four values of  in the mantle and n = 3. For a given Newtonian viscosity, non-linear rheology dominates at a certain location if the von Mises stress is locally higher than that given by the curve.

until the end of deglaciation—this is generally supported by the findings of Wu (1999) and Wu and Wang (2008). Farther away from the ice center (e.g. around 10.25◦ ), non-linear creep is still smaller than linear creep in the upper mantle. However, near the ice margin (e.g. around 15.25◦ ), q >qt between 18 and 8 ka BP, thus a transition from non-linear to linear creep as the dominant mechanism occurs around 8 ka BP. 4. Misfit and relative sea level curves

Fig. 2. Time evolution of von Mises stress q at four locations just below the lithosphere (at r = 6220 km) in the axisymmetric model of Fig. 1. The horizontal thick line indicates the transitional von Mises stress level (see Fig. 3), above which the dominant creep is non-linear.

Model predictions are compared with RSL data for 128 sites that are taken from the Tushingham and Peltier (1991) database. These 128 sites are selected with the following criteria: sites with less than 4 data points and spanning less than 4 ka and sites which cannot clearly indicate a consistent trend were removed. In addition, inconsistent and problematic sea level data points were also removed from the remaining sites. We have tested the sensitivity of these 128 sites by removing the site that has the largest misfit among all the sites and verified that the best-fitting model (linear versus non-linear versus composite) and the best-fitting model

42

W. van der Wal et al. / Journal of Geodynamics 50 (2010) 38–48

Fig. 4. Misfit between 128 RSL curves and predictions with the ICE-5G loading history for linear, non-linear and composite rheologies for four values of the Newtonian viscosity: (a) 1 × 1021 Pa s, (b) 3 × 1021 Pa s, (c) 9 × 1021 Pa s and (d) 27 × 1021 Pa s.

parameters remained the same. All ages are C14 age calibrated according to Fairbanks et al. (2005). We define misfit:

  1  oi − pi 2

F=

N

2

 =

N

i=1

oi

statistically significant we have used the F-test, following Dal Forno and Gasperini (2007). The F-statistic is computed as

(4)

where N is the number of observations (1071 RSL observations), oi are the observations, pi are the predicted values from the models, and  oi are the standard deviations of the observations (being a combination of time and height error expressed in RSL height). The misfit for linear, non-linear and composite models is plotted in Fig. 4a for four different values of A* and a Newtonian viscosity of 1 × 1021 Pa s in combination with the ICE-5G loading history. In Fig. 4b, c and d the Newtonian viscosity is increased to 3 × 1021 , 9 × 1021 and 27 × 1021 Pa s, respectively. Except for the case of  = 1 × 1021 , the best-fitting rheology is always the composite rheology. In order to check if the improvement in misfit is

(21 − 22 )/( 2 − 1 ) 22 /(N − 2 )

,

(5)

where subscript 2 denotes the composite rheology, subscript 1 denotes the Chi-squared of the linear or non-linear rheology, are the degrees of freedom in the model (3 for composite, 2 for nonlinear and 1 for linear rheology) and N is the number of observations (1071). The improved fit of the best-fitting composite rheology for the Newtonian viscosities in Fig. 4b, c and d is statistically significant at the 99% level compared to the fit of the linear and non-linear rheology models. The pre-exponent factor that provides the best fit is A* = 3.3 × 10−35 Pa−3 s−1 in all cases (as in Wu, 1999). For later reference, we label the best-fitting composite rheology model with A* = 3.3 × 10−35 Pa−3 s−1 and  = 3 × 1021 Pa s, ‘Model I’, and the second best-fitting model, with A* = 3.3 × 10−35 Pa−3 s−1 and

Fig. 5. Location of the 12 RSL sites used for Fig. 6.

W. van der Wal et al. / Journal of Geodynamics 50 (2010) 38–48

43

Fig. 6. RSL predictions from models with the ICE-5G loading history and with three different rheologies: linear, non-linear and composite rheology, for the 12 RSL sites of Fig. 5. The stress exponent n = 3, A* = 3.3 × 10−35 Pa−3 s−1 and  = 3 × 1021 Pa s.

 = 9 × 1021 Pa s, ‘Model II’. Transition stresses for these models are 1.8 and 0.6 MPa, respectively. The misfits of Model I (6.5) and Model II (7.0) are close. A Newtonian viscosity of 3 × 1021 Pa s is close to the Newtonian viscosity of 2.7 × 1021 Pa s for the best-fitting composite rheology in Dal Forno and Gasperini (2007). With a transition stress of 1.5 MPa the pre-exponent factor A* for their best-fitting model becomes 1.8 × 10−35 Pa−3 s−1 , which is also close to Model I. Finally, note that the misfit curves for the composite rheology are generally close to the curves for non-linear rheology. This is due to the low transition stress of Model I which makes non-linear rhe-

ology dominant in and around the centers of rebound most of the time. Misfit provides only one number for the entire set of RSL sites, while sometimes the detailed RSL response at the ice margin is useful as distinction between linear and non-linear rheology (Wu, 1995). Therefore, in the following we compare the RSL curves predicted by the composite rheology with those predicted by purely non-linear and purely linear rheologies, for a selection of 12 sites depicted in Fig. 5. RSL curves are plotted in Fig. 6 for Model I (A* = 3.3 × 10−35 Pa−3 s−1 ,  = 3 × 1021 Pa s) as well as linear and non-linear rheology models with the same parameters. The transition stress is low enough for non-linear deformation to be dominant near the centers of rebound so that the RSL curve for composite rheology is mostly close to the curve for the purely non-linear case. However, at sites near the margin of the Laurentide ice sheet (Newfoundland, Boston—up to 12 ka BP) and in the center and margin of the Fennoscandian ice sheet (Bjugn, Onsala, Angermanland, Helsinki and Kong Karlsland) RSL curves of the composite rheology drop faster than non-linear curves, apparently because of the transition to linear deformation as shown in Fig. 2. Note that for Newfoundland this contribution leads to a better fit than either linear or non-linear rheology. An exception is Moruya, where the stress is apparently not large enough to trigger large non-linear deformation. 5. Uplift rate

Fig. 7. Maximum uplift rate in North America for non-linear and composite rheologies for different Newtonian viscosities, different values of A* with the ICE-5G loading history.

The best-fitting model to RSL curves among the limited number of models investigated in this paper is that with composite rheology. However, a problem with non-linear rheologies is the low predicted present-day uplift rates (e.g. Wu, 1999). Although the relative sea level curves are the most important for GIA inferences because their records extend back in time, an uplift rate that is too low points to a problem in the model or its inputs. It is useful to see

44

W. van der Wal et al. / Journal of Geodynamics 50 (2010) 38–48

Fig. 8. (a) Maximum uplift rate in North America for four composite rheologies with (delayed versions of) ICE-4G. Solid line denotes no delay; dashed line denotes a delay of 2 ka. (b–d) RSL misfit with respect to those sites that are located on the North American continent for composite rheologies with (delayed versions of) ICE-4G, with Newtonian viscosity of 1 × 1021 Pa s (b), 3 × 1021 Pa s (c) and 27 × 1021 Pa s (d).

if composite rheology models might be able to provide larger uplift rates in addition to an improved fit with RSL data. Following Wu and Wang (2008), we focus on the maximum uplift rate (in North America) because the uplift rate pattern is determined more by the ice model. As reviewed in Wu and Wang (2008), the maximum observed uplift rate in Laurentide is 11 ± 2 mm/a, and this is to be compared with the predicted maximum uplift rates in Fig. 7. Note that ICE-5G is used here. Inspection of the figure shows that indeed composite rheology can result in a higher uplift rate than that predicted by a purely non-linear rheology, and the predicted rate is highest when Newtonian viscosity is around 3 (Model I) to 9 × 1021 Pa s (Model II). Maximum uplift rate of 7.5 mm/a is found for Model I and Model II. This uplift rate is far too small, considering the 11 mm/a cited above and considering that the ICE-5G model already gives larger uplift rates than the ICE-4G model. However, both global ice histories ICE-4G and ICE-5G are constructed with linear rheology (Peltier, 1994, 2004), therefore these histories should be modified if mantle rheology is non-linear or composite. 6. Modification of ice history In this section, we investigate briefly the modification of ice history required to better fit the observed uplift rate. Because of the long computation time of FE models with ABAQUS we are limited to investigating a few simple cases of modification to the ice history in North America. In the following, we investigate the effect of a time delay in glaciation (which has been shown to produce higher uplift rates in Wu and Wang, 2008) and the effect of an increase in ice height in North America (as suggested by Wu and Wang, 2008), on GIA observables with composite rheology. These modifications merely serve as directions for producing future ice models based on non-linear or composite rheology, we do not claim that the delay and scaling of the ice height are allowed by geological constraints

or can fit GIA constraints other than the ones we consider. Glaciologists found that the ice thickness in ICE-4G is ∼1.5 times too thin (Hughes, 1998) and this led Peltier to develop the ICE-5G model. Unfortunately, the predicted uplift rates from ICE-5G are somewhat worse than that of ICE-4G (van der Wal et al., 2009). We have chosen to use the older ICE-4G model (Peltier, 1994) to demonstrate the effect of rheology on uplift rates, even though ICE-5G is in many aspects a significant improvement to the ICE-4G model. In the following, maximum uplift rates will be discussed, as well as misfit values and sea level curves. First we show the maximum uplift rate achieved by different composite rheology models in Fig. 8a. Uplift rates increase with decreasing A* for all values of the Newtonian viscosity. A delay in glaciation increases the uplift rate for all models, but slightly less so for models with a larger value of A* which have a larger nonlinear component and which tend to have flatter sea level curves around present. The misfit values for RSL data in North America are plotted in Fig. 8b, c and d for three different values of the viscosity. Again it was verified that removing the site with the largest misfit did not significantly change the results. The difference between  = 3 × 1021 and  = 9 × 1021 Pa s is small, therefore the latter is not shown. A delay in the ICE-4G glaciation history leads to an increase in misfit of RSL data for small values of A* and less bad or even better fit for A* = 3.3 × 10−34 Pa−3 s−1 . The improved misfit for 1 ka delay for A* = 3.3 × 10−34 Pa−3 s−1 and for all viscosities is statistically significant at the 99% level. The best-fitting model is found for no delay, A* = 3.3 × 10−35 Pa−3 s−1 and  = 3 × 1021 Pa s (Model I). However, a 1 ka delay seems to be acceptable based on the sea level data. With a 1 ka delay in glaciation we can increase uplift rate for this model from 5.3 to approximately 6.4 mm/year (halfway between the values for 0 and 2 ka delay) while the RSL misfit only increases from 4.7 to 5.0. To have a visual comparison of the fit with relative sea level data, Fig. 9 plots the RSL curves at 12 sites in North America (see

W. van der Wal et al. / Journal of Geodynamics 50 (2010) 38–48

45

Fig. 9. RSL predictions for the 12 North American sites in Fig. 10, with delayed versions of ICE-4G and for a composite rheology with n = 3, A* = 3.3 × 10−35 Pa−3 s−1 and  = 3 × 1021 Pa s.

Fig. 10) for Model I with 0, 1 and 2 ka delay applied to ICE-4G. It can be seen that a 2 ka delay improves the fit for many sites with land emergence (Richmond Gulf, Churchill, Southampton Island, Milne Inlet Baffin Island, Bathurst Inlet Northwest Territories), confirming the conclusion of Wu and Wang (2008) for non-linear rheology. A smaller delay of 1 ka seems to better fit sites like Ottawa Island, North Somerset Island, Newfoundland, and to a lesser extent Clinton and Brigantine. To investigate the effect of scaling ice heights on uplift rate and RSL predictions we scale the ICE-4G ice heights in Laurentide (as

suggested by Wu and Wang, 2008 and also done by Wu, 1999) by multiplying ice heights at all time steps by a factor of 1.5 and 2.0. In doing so we violate the far-field constraints on total ice volume, but here we are only interested in whether a change in ice thickness in combination with composite rheology can provide a better fit to uplift rate data and sea level data in North America. Fig. 11 shows the maximum uplift rate achieved with these scaled versions of ICE-4G. Doubling the ice heights adds 1.8 mm/a to the maximum uplift rate for Model I and 1.7 mm/a for Model II. Because of the lower transition stress, non-linear creep dominates and relaxation proceeds slightly faster for Model II—thus the effect of increasing ice thickness on uplift rate is slightly smaller. Fig. 11b shows a worse sea level fit for all values of A*. Misfit increases less when the viscosity is increased to 3 and 27 × 1021 Pa s. For Model II scaling has little effect on the misfit, but uplift rate is also not increased significantly. One solution to increase uplift rate at the expense of a slightly lower RSL misfit value is to decrease the value of A* compared to Model I (see also Fig. 8) and perhaps increase  to a value in between 3 and 9 × 1021 Pa s. In Fig. 12 the sea level curves for Model I with increased ice thickness are shown to be closer to observations at some sites within and close to the ice margin (Ottawa Island, Milne Inlet, North Somerset Island, Boston). This agrees with the finding of Wu (1999). However, the shapes of the sea level curves for upscaled ice heights in sites Richmond Gulf, Churchill and Southampton Island do not fit the observations as well as the sea level curves resulting from a delay in glaciation do. 7. Discussion

Fig. 10. Location of the North American RSL sites used for Figs. 9 and 12.

A composite rheology in which both diffusion and dislocation creep operate at the same time is more in line with microphysics of mantle creep. Here we found that GIA models based on composite rheology with mantle averaged values of A* and  show a

46

W. van der Wal et al. / Journal of Geodynamics 50 (2010) 38–48

Fig. 11. Maximum uplift rate in North America for four composite rheologies with (scaled versions) of ICE-4G. Solid line denotes no scaling, dashed line denotes scaling by 1.5, dotted line denotes scaling by 2.0. (b–d) RSL misfit for composite rheologies with (scaled versions of) ICE-4G, with Newtonian viscosity of 1 × 1021 Pa s (b), 3 × 1021 Pa s (c) and 27 × 1021 Pa s (d).

Fig. 12. RSL predictions for the 12 North American sites in Fig. 10, using the ICE-4G model with increased ice heights. Composite rheology is used with n = 3, A* = 3.3 × 10−35 Pa−3 s−1 and  = 3 × 1021 Pa s.

W. van der Wal et al. / Journal of Geodynamics 50 (2010) 38–48

statistically significant better fit to RSL data than models based on purely linear or non-linear rheology (in agreement with Gasperini et al., 2004) for the 16 models investigated here. The thickness of the lithosphere is an important issue, especially for the shape of the peripheral bulge, which is not discussed in this study. At first sight, it seems as if composite rheology has the advantage of an extra parameter compared to linear or non-linear rheology and should therefore always show the smallest misfit among simulations performed over the entire parameter space. However, this is not necessarily the case, as the relative importance of both parameters is coupled through the state of stress. As an example, consider the best-fitting non-linear rheology model (A* = 3.3 × 10−35 Pa−3 s−1 ;  = 3 × 1021 Pa s) and the composite rheology models with  = 1 × 1021 Pa s. The latter have sea level curves closer to those of linear rheology which gives higher RSL misfit and is worse than purely non-linear rheology for the values of A* investigated here as can be seen in Fig. 4a. The conclusions regarding the research questions posed in the Introduction are presented below: (1) It has been shown that the region where non-linear creep dominates over linear creep occurs in an area beneath and close to the former ice sheet margin which grows in size until the end of melting. (2) The best-fitting composite rheology models behave mostly in a non-linear way but at locations near the Laurentide ice margin and the center and edge of the Fennoscandian ice sheet the contribution from linear rheology becomes more important. In view of this, results from previous studies with a purely linear or purely non-linear rheology in (parts of) the mantle are still relevant even if the real Earth responds to glacial loading with a composite rheology, because the best-fitting composite rheology behaves as a linear or non-linear rheology for part of the deglaciation. (3) A composite rheology with a Newtonian viscosity of 3 × 1021 Pa s and pre-stress exponent of 3.3 × 10−35 Pa−3 s−1 Pa s has the lowest RSL misfit of all models investigated here (for the ICE-4G and ICE-5G ice loading histories), and at the same time slightly increases the present-day uplift rate compared to a purely non-linear rheology. To a lesser extent this also holds for a Newtonian viscosity of 9 × 1021 Pa s. This is an encouraging result that supports the use of composite rheology in future GIA modeling. A larger uplift rate can probably also be obtained for values of the Newtonian viscosity between 3 and 9 × 1021 Pa s and a value of the pre-stress exponent between 3.3 × 10−35 and 3.3 × 10−36 Pa−3 s−1 , at the expense of a slightly larger RS misfit. (4) Both increasing ice thickness and delaying the glacial cycle lead to higher uplift rates, but less so for increased values of A*. Delay in glaciation can lead to a decrease in RSL misfit which is statistically significant, but only for a value of A* of ∼3.3 × 10−34 Pa−3 s−1 where the uplift rate is far too small. Both scaling and delay improve fit with sea level curves at some sites but perform worse in terms of overall misfit. Delay in glaciation can be expected to be closer to field observations than increase in ice height, while scaling the ice heights in North America by 1.5 would require much thinner ice in the rest of the world in order to fit far-field constraints on the total melt-water volume. Thus, delay in glaciation seems a more promising direction for tailoring an ice model to composite rheology. The importance of linear and non-linear deformation mechanisms depends strongly on many parameters such as activation volume, activation enthalpy, grain size and water content (Karato and Wu, 1993) which are grouped in the pre-stress exponent A* which is assumed to be constant throughout the mantle in this

47

paper. Composite rheology can be implemented in multi-layer Earth models in the future. However, it remains to be seen if GIA data has enough resolving power to justify inference of creep parameters in multi-layer Earth. Acknowledgements The authors would like to thank Doug Philips for computer help, Paolo Gasperini and Thorsten Becker for comments on an earlier version of this manuscript and Volker Klemann, an anonymous reviewer and guest editor Markku Poutanen for helping to improve this manuscript. Patrick Wu is supported by a Discovery Grant from NSERC (Canada). Hansheng Wang is supported by the National Natural Science Fund of China for Distinguished Young Scholars (40825012). Michael Sideris is supported by the GEOIDE Network of Centers of Excellence (Canada). References Bird, R.B., Stewart, W.E., Lightfoot, E.N., 1960. Transport Phenomena. Wiley, New York. Cordier, P., Ungar, T., Zsoldos, L., Tichy, G., 2004. Dislocation creep in MgSiO3 perovskite at conditions of the Earth’s uppermost lower mantle. Nature 428, 837–840. Dal Forno, G., Gasperini, P., Boschi, E., 2005. Linear or nonlinear rheology in the mantle: a 3D finite-element approach to postglacial rebound modeling. J. Geodyn. 39, 183–195. Dal Forno, G., Gasperini, P., 2007. Modeling of mantle postglacial relaxation in axisymmetric geometry with a composite rheology and a glacial load interpolated by adjusted spherical harmonics analysis. Geophys. J. Int. 169, 1301–1311, doi:10.1111/j.1365-246X.2007.03347.x. Fairbanks, R.G., Mortlock, R.A., Chiu, T.C., Cao, L., Kaplan, A., Guilderson, T.P., Fairbanks, T.W., Bloom, A.L., Grootes, P.M., Nadeau, M.-J., 2005. Radiocarbon calibration curve spanning 0 to 50,000 years BP based on paired 230 Th/234 U/238 U and 14 C dates on pristine corals. Quat. Sci. Rev. 24, 1781–1796, doi:10.1016/j.quascirev.2005.04.007. Gasperini, P., Yuen, D.A., Sabadini, R., 1992. Postglacial rebound with a nonNewtonian upper mantle and a Newtonian lower mantle rheology. Geophys. Res. Lett. 19, 1711–1714. Gasperini, P., Dal Forno, G., Boschi, E., 2004. Linear or nonlinear rheology in the Earth’s mantle: the prevalence of power-law creep in the postglacial isostatic readjustment of Laurentia. Geophys. J. Int. 157, 1297–1302. Giunchi, C., Spada, G., 2000. Postglacial rebound in a non-Newtonian spherical Earth. Geophys. Res. Lett. 27, 2065–2068. Goetze, C., Kohlstedt, D.L., 1973. Laboratory study of dislocation climb and diffusion in olivine. J. Geophys. Res. 78, 5961–5971. Hughes, T., 1998. Tutorial on strategies for using isostatic adjustments in models that reconstruct ice sheets during the last deglaciation. In: Wu, P. (Ed.), Dynamics of the Ice Age Earth: A Modern Perspective,. Trans Tech Publ., Switzerland, pp. 365–382. Karato, S.-I., Wu, P., 1993. Rheology of the upper mantle: a synthesis. Science 260, 771–778. Karato, S.-I., 1998. Micro-physics of post-glacial rebound. In: Wu, P. (Ed.), Dynamics of the Ice Age Earth: A Modern Perspective. Trans Tech Publ., Switzerland, pp. 351–364. Korenaga, J., Karato, S.-I., 2008. A new analysis of experimental data on olivine rheology. J. Geophys. Res. 113, B02403, doi:10.1029/2007JB005100. Li, L., Raterron, P., Weidner, D., Chen, J., 2003. Olivine flow mechanisms at 8 GPa. Phys. Earth Planet. Int. 138, 113–129. Mainprice, D., Tommasi, A., Couvy, H., Cordier, P., Frost, D.J., 2005. Pressure sensitivity of olivine slip systems and seismic anisotropy of Earth’s upper mantle. Nature 433, 731–733. Parmentier, E.M., Turcotte, D.L., Torrance, K.E., 1976. Studies of finite amplitude nonNewtonian thermal convection with application to convection in the Earth’s mantle. J. Geophys. Res. 81, 1839–1846. Peltier, W.R., 1984. The thickness of the continental lithosphere. J. Geophys. Res. 89, 11,303–11,316. Peltier, W.R., 1994. Ice age paleotopography. Science 265, 195–201. Peltier, W.R., 1998. Postglacial variations in the level of the sea: implications for climate dynamics and solid-Earth geophysics. Rev. Geophys. 36, 603– 689. Peltier, W.R., 2004. Global glacial isostasy and the surface of the ice-age earth: the ICE-5G (VM2) Model and Grace. Annu. Rev. Earth Planet. Sci. 32, 111– 149. Ranalli, G., 1991. The microphysical approach to mantle rheology. In: Sabadini, R., Lambeck, K., Boschi, E. (Eds.), Glacial Isostacy, Sea level and Mantle Rheology. Kluwer Academic Publishers, Dordrecht, pp. 343–378. Ranalli, G., 1995. Rheology of the Earth, 2nd ed. Chapman and Hall, London. Ranalli, G., 2001. Mantle rheology: radial and lateral viscosity variations inferred from micro-physical creep laws. J. Geodyn. 32, 425–444.

48

W. van der Wal et al. / Journal of Geodynamics 50 (2010) 38–48

Schmeling, H., 1987. On the interaction between small- and large-scale convection and postglacial rebound flow in a power-law mantle. Earth Planet. Sci. Lett. 84, 254–262. Turcotte, D.L., Schubert, G., 2002. Geodynamics, 2nd ed. Cambridge University Press, Cambridge. Tushingham, A., Peltier, W.R., 1991. ICE-3G: a new global model of late Pleistocene deglaciation based upon geophysical predictions of post-glacial relative sea level change. J. Geophys. Res. 96, 4497–4523. van den Berg, A., van Keken, P.E., Yuen, D.A., 1993. The effects of a composite nonNewtonian and Newtonian rheology on mantle convection. Geophys. J. Int. 115, 62–78. Van der Wal, W., Braun, A., Wu, P., Sideris, M.G., 2009. Prediction of decadal slope changes by glacial isostatic adjustment modeling. Can. J. Earth Sci. 46, 587–595. Wang, H., Wu, P., 2006. Effects of lateral variations in lithospheric thickness and mantle viscosity on glacially induced surface motion on a spherical, self-gravitating Maxwell Earth, Earth Planet. Sci. Lett. 244, 576– 589. Wolf, D., Klemann, V., Wunsch, J., Zhang, F.-P., 2006. A reanalysis and reinterpretation of geodetic and geological evidence of glacial-isostatic adjustment in the Churchill region, Hudson Bay. Surv. Geophys. 27, 19–61, doi:10.1007/s10712005-0641-x.

Wu, P., 1995. Can observations of postglacial rebound tell whether the rheology of the mantle is linear or nonlinear? Geophys. Res. Lett. 22, 1645–1648. Wu, P., 1998. Postglacial rebound modeling with power law rheology. In: Wu, P. (Ed.), Dynamics of the Ice Age Earth: A Modern Perspective. Trans Tech Publ., Switzerland, pp. 365–382. Wu, P., 1999. Modelling postglacial sea levels with power-law rheology and a realistic ice model in the absence of ambient tectonic stress. Geophys. J. Int. 139, 691–702. Wu, P., 2001. Postglacial induced surface motion and gravity in Laurentia for uniform mantle with power-law rheology and ambient tectonic stress. Earth Planet. Sci. Lett. 186, 427–435. Wu, P., 2002. Effects of mantle flow law stress exponent on postglacial induced surface motion and gravity in Laurentia. Geophys. J. Int. 148, 676–686. Wu, P., van der Wal, W., 2003. Postglacial sea levels on a spherical, self-gravitating viscoelastic earth: effects of lateral viscosity variations in the upper mantle on the inference of viscosity contrasts in the lower mantle. Earth Planet. Sci. Lett. 211, 57–68. Wu, P., 2004. Using commercial finite element packages for the study of Earth deformations, sea levels and the state of stress. Geophys. J. Int. 158, 401–408. Wu, P., Wang, H., 2008. Postglacial isostatic adjustment in a self-gravitating spherical Earth with power-law rheology. J. Geodyn. 46, 118–130, doi:10.1016/ j.jog.2008.03.008.