# Effects of mantle rheologies on viscous heating induced by Glacial Isostatic Adjustment

Effects of mantle rheologies on viscous heating induced by Glacial Isostatic Adjustment Summary It has been argued that viscous dissipation from mantle flow in response to surface loading during glacial cycles can result in short-term heating and thus trigger transient volcanism or changes in mantle properties, which may in turn affect mantle dynamics. Furthermore, heating near the Earth's surface can also affect the stability of ice sheets. We have studied the magnitude and spatial-temporal distribution of viscous heating induced in the mantle by the realistic ice model ICE-6G and gravitationally consistent ocean loads. Three types of mantle rheologies, including linear, non-linear and composite rheologies are considered to see if non-linear creep can induce larger viscous heating than linear rheology. We used the Coupled-Laplace-Finite-Element model of Glacial Isostatic Adjustment (GIA) to compute the strain, stress and shear heating during a glacial cycle. We also investigated the upper bound of temperature change and surface heat flux change due to viscous heating. We found that maximum viscous heating occurs near the end of deglaciation near the edge of the ice sheet with amplitude as high as 120 times larger than that of the chondritic radioactive heating. The maximum heat flux due to viscous heating can reach 30 mW m−2, but the area with large heat flux is small and the timescale of heating is short. As a result, the upper bound of temperature change due to viscous heating is small. Even if 30 glacial cycles are included, the largest temperature change can be of the order of 0.3 °C. Thus, viscous heating induced by GIA cannot induce volcanism and cannot significantly affect mantle material properties, mantle dynamics nor ice-sheet stability. Loading of the Earth, Europe, North America, Heat generation and transport, Rheology: mantle 1 INTRODUCTION In a fluid, the presence of viscosity, which represents the frictional interaction between neighbouring fluid elements, transforms some kinetic energy of the fluid into heat energy for dissipation (e.g. Landau & Lifshitz 1966). This so-called ‘viscous heating’ can potentially change the temperature and therefore material properties of mantle rocks. If viscous heating is large, then it can be very important because it can induce melting and therefore volcanism. It can also alter seismic properties and thus our interpretation of internal structure of the Earth, or alter mantle viscosity, which may in turn affect mantle dynamics. For example, viscous heating generated by tidal deformation induced by Jupiter may have caused volcanism on the Jovian moon Io (Ross & Schubert 1987; Segatz et al.1988). There is speculation that glacial cycles may have triggered volcanism during the Quaternary (e.g. Nakada & Yokose 1992; Sigvaldason et al.1992; Huybers & Langmuir 2009, Uenzelmann-Neben et al.2012). The link between them is generally attributed to decompression (Jull & McKenzie 1996; Slater et al.1998), but volcanism can also be triggered by viscous heating induced by Glacial Isostatic Adjustment (GIA) of the Earth in response to cycles of ice and ocean loading. Hanyk et al. (2005) modeled viscous heating due to GIA using simple ice load histories on a spherical Earth with linear rheology. They found that viscous heating can generate observable transient heat flux on the surface of the Earth. They also suggested that some degree of volcanism may be associated with viscous heating if the deglaciation period is short, and speculated that higher values of viscous heating can be generated if mantle rheology is non-linear. Indeed, if viscous heating is large enough at the base of the ice sheet, then glacial loading can also affect the stability of ice sheets by basal heating (e.g. Pattyn 2010). On the other hand, if viscous heating leads to lower mantle viscosity near the surface, then this would lead to faster land uplift and that may help to stabilize marine ice sheets (Gomez et al.2015). In this paper, we improve on Hanyk's work by using a realistic ice history model ICE-6G (Peltier et al.2015) together with realistic, self-gravitating oceans. As we shall see, the interaction among ice sheets and the loading of the ocean floor by meltwater affects the spatial distribution of viscous heating. A realistic loading history is also important, as it will eliminate the unrealistic results due to unreasonably fast or slow deglaciation histories. In addition, we consider three types of mantle rheologies in our Earth models: linear, non-linear and composite rheologies to see how large non-linear creep can affect the magnitude of viscous dissipation. Besides studying the effect of viscous heating on surface heat flux, we also study the following questions: what temperature changes inside the Earth are produced by viscous heating induced by glacial cycles? Do these temperature changes significantly affect present-day seismic velocities and thus the interpretation of seismic tomography? Can they affect viscosity and thus mantle flow induced during GIA? In the following sections, we will start with a brief review of the three different types of mantle rheologies, and provide a brief discussion of our model. Then, we will present the results, focusing on the magnitude and spatial and temporal distribution of viscous heating for some simple and other more sophisticated Earth models. After that, we will estimate the effects of viscous heating on the upper bound of heat flux change and see if the upper bound of temperature increase can trigger volcanism and affect ice stability or material properties inside the Earth. 2 THE MODEL 2.1 Mantle rheologies and rate of viscous dissipation Creep experiments on mantle rocks show that both linear and non-linear creep laws operate in the mantle (e.g. Karato & Wu 1993). If both creep laws operate simultaneously, then we have composite rheology. For this case, the creep mechanism with the highest creep rate becomes the dominant creep mechanism. The constitutive relation for these rheologies relates the deviatoric stress σD and the deviatoric strain rate $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$, which are defined as   $${{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}} = {\boldsymbol{\sigma }} - \frac{1}{3}{\sigma _{ii}}{{\bf I}}$$ (1)  $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}} = {\boldsymbol{\dot{\varepsilon }}} - \frac{1}{3}{\dot{\varepsilon }_{ii}}{{\bf I}}$$ (2)where σ and $${\boldsymbol{\dot{\varepsilon }}}$$ are the stress tensor and the strain rate tensor, respectively, I is the identity tensor, σii is the sum of diagonal components of σ and $${\dot{\varepsilon }_{ii}}$$ is the sum of diagonal components of $${\boldsymbol{\dot{\varepsilon }}}$$. In the following, we refer to σDij and $${\dot{\varepsilon }_{Dij}}$$ as the components of σD and $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$, respectively. The relation between σD and $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ for composite, linear and non-linear rheologies can be expressed as (e.g. van der Wal et al.2010):   $${{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}} = \frac{1}{{\frac{1}{{2\eta }} + A{\sigma _E}^{n - 1}}}{{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ (3a) Here, η is the dynamic viscosity for linear creep, A is the non-linear creep parameter determined from shear experiments, and its value depends on temperature, pressure and material properties. σE is the effective deviatoric stress defined by $${\sigma _E} = \sqrt {\frac{1}{2}{\sigma _{Dij}}{\sigma _{Dij}}}$$ and n is the stress exponent with experimental value between 2 and 6. Here, we will take n = 3, as that is a typical value for mantle rocks (Karato & Wu 1993, Ranalli 1995). If A = 0, then eq. (3a) reduces to the relation for linear rheology. If A ≠ 0 and η = ∞, then it becomes the relation for non-linear rheology. It is useful to define the effective viscosity $${\eta _{{\rm{eff}}}} = \frac{1}{2}{( {\frac{1}{{2\eta }} + A{\sigma _E}^{n - 1}} )^{ - 1}}$$, so that eq. (3a) becomes   $${{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}} = 2{\eta _{{\rm{eff}}}}{{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ (3b) For non-linear or composite rheology, the effective viscosity changes with stress level and thus with space and time. For non-linear rheology or composite rheology with a fixed η, a decrease in σE will result in an increase in effective viscosity. The viscous dissipation rate ϕ (hereafter called ‘viscous heating’) is given by (e.g. eq. 7 in Hanyk et al.2005):   $$\phi = {{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}}{:}\, {{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}} = \frac{1}{{2{{\boldsymbol{\eta }}_{{{\bf eff}}}}}}{{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}}{:}\,{{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}}$$ (4a) This shows that shear heating decreases if the effective viscosity increases provided that σD does not change rapidly. For surface loading problems like GIA, σD is determined by the changing surface loads, although σD is also affected by stress relaxation where the decay time constant is proportional to the effective viscosity. In terms of von Mises stress $$\tau = \sqrt 3 {\sigma _E}$$, the creep parameter A and viscosity η, eq. (4a) can be written as:   $$\phi = \frac{2}{3} \left( {\frac{1}{{2\eta }} + \frac{1}{3}A{\tau ^2}} \right){\tau ^2}$$ (4b)which clearly shows the strong dependence of ϕ on von Mises stress τ. Eq. (4a) can also be written as: $$\phi = 2{{\boldsymbol{\eta }}_{{{\bf eff}}}}{{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}{:}\,{{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ which shows that viscous heating would increase with larger effective viscosity provided that the strain rate $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ is constant. For our problem, $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ is not constant, but $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}{:}\,{{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ is determined by (σD/ηeff)2. Thus, eqs (4a) or (4b) will be used for our discussion below. 2.2 The GIA model The Coupled Laplace-Finite-Element (CLFE) GIA model of Wu (2004), later modified by van der Wal et al. (2010), is used in this study. Finite-element grids with different (0.5° and 2°) spatial resolution have been used for the computations. It is found that 2° resolution is adequate for our purpose. The inputs of the GIA model are elastic and creep parameters for the Earth model, and ice loading history. Unlike the simple ice sheet model used in Hanyk et al. (2005), the ice model used here is the realistic global ICE6G model of Peltier et al. (2015). Since ICE6G provides ice thickness history from 26 kBP to the present, we assume that the ice thickness increased linearly from zero at 108 kBP to the ice thickness at 26 kBP. Also, the ice was taken from the water in the oceans and meltwater returned to the oceans. The self-consistent sea level equation is solved for realistic oceans (Wu 2004). Since the effects of time-dependent coastline and rotational feedback are small on the von Mises stress τ and shear heating, they have not been included. As we shall see below, the localization of viscous heating in time makes the consideration of previous glacial cycles unnecessary. For the Earth models, the elastic parameters for models M1–M3 are the same as those in van der Wal et al. (2010), and for M4, they are the same as that in model VM5a of Peltier et al. (2015). The rheological parameters are described in Table 1. Here, parameter A* is the creep parameter determined from uniaxial experiments and is related to A, the creep parameter from shear experiments by $${A^*} = 2A/\sqrt {{3^{n + 1}}}$$ (van der Wal et al.2010). In Table 1, M1 and M4 are linear rheological models. M1 has a uniform viscosity of 3 × 1021 Pa s, while in M4, the viscosities in the upper mantle (100–670 km depth), shallow lower mantle (670–1271 km depth) and deep lower mantle (1271–2891 km depth) are 5 × 1020, 1.6 × 1021 and 3 × 1021 Pa s respectively. M4 is considered because model VM5a is used for the construction of ICE6G (Peltier et al.2015), and so these should be used together. The rheological parameters in models M2 and M3 are uniform in the mantle. M2 has a non-linear rheology, while model M3 has a composite rheology. In other models, the effects of compressibility and a low viscosity (LV) zone below the lithosphere are also studied (see Tables 2–4) as they can also affect the magnitude of viscous heating. Table 1. Rheological parameters of four Earth models.   M1  M2  M3  M4  A*(Pa − 3 · s − 1)  0  1.11 × 10−34  1.11 × 10−34  0  η (Pa s)  3 × 1021  0  3 × 1021  VM5a  n  1  3  3  1    M1  M2  M3  M4  A*(Pa − 3 · s − 1)  0  1.11 × 10−34  1.11 × 10−34  0  η (Pa s)  3 × 1021  0  3 × 1021  VM5a  n  1  3  3  1  View Large Table 2. The maximum local viscous heating of all time for uniform mantle models with various rheology (A* and η). ϕMax  η(Pa · s)      3.00 × 1020  3.00 × 1021  3.00 × 1022  Non-linear  A* (Pa − 3 · s − 1)  Linear  11.64  3.95  0.54      1.11 × 10−36  11.64  6.14  2.73  2.23    1.11 × 10−35  11.45  11.45  10.24  9.99    1.11 × 10−34  9.58  10.04  10.12  10.14    1.11 × 10−33  6.73  6.55  6.54  6.53  ϕMax  η(Pa · s)      3.00 × 1020  3.00 × 1021  3.00 × 1022  Non-linear  A* (Pa − 3 · s − 1)  Linear  11.64  3.95  0.54      1.11 × 10−36  11.64  6.14  2.73  2.23    1.11 × 10−35  11.45  11.45  10.24  9.99    1.11 × 10−34  9.58  10.04  10.12  10.14    1.11 × 10−33  6.73  6.55  6.54  6.53  View Large Table 3. The maximum local viscous heating of all epochs for two values of Poisson's ratio (i.e. compressibility) with various rheology (A* [Pa−3 s−1] and η [Pa s]). ϕMax  A* = 0 η = 3 × 1021  A* = 1.11 × 10−34 (non − linear)  A* = 1.11 × 10−34η = 3 × 1021  Poisson's ratio  0.4900  3.95  10.14  10.04    0.2877  6.18  10.59  10.40  ϕMax  A* = 0 η = 3 × 1021  A* = 1.11 × 10−34 (non − linear)  A* = 1.11 × 10−34η = 3 × 1021  Poisson's ratio  0.4900  3.95  10.14  10.04    0.2877  6.18  10.59  10.40  View Large Table 4. The maximum local viscous heat of all time in models with an LV zone below the lithosphere. Units: A∗ in [Pa−3 s−1] and η in [Pa s]. ϕMax  A*=0 η = 3 × 1018  A* = 1.11 × 10−34 (non − linear)  A* = 1.11 × 10−36η = 3 × 1018  LV zone thickness (km)  40  102.18  26.01  102.19    100  118.80  30.52  118.81  ϕMax  A*=0 η = 3 × 1018  A* = 1.11 × 10−34 (non − linear)  A* = 1.11 × 10−36η = 3 × 1018  LV zone thickness (km)  40  102.18  26.01  102.19    100  118.80  30.52  118.81  View Large The outputs of the CLFE model contain spatial temporal evolution of the displacements and state of stress throughout the mantle. From the stress output at any time step, the local viscous heating for each element can be computed using eqs (4a) or (4b). In this paper, viscous heating is normalized by the chondritic radiogenic heating of 3 × 10−9 W m−3 (Hanyk et al.2005). 3 RESULTS 3.1 Spatial-temporal distribution of viscous heating for M1–M4 In Fig. 1, the spatial distribution of normalized viscous heating for models M1–M4 is shown at 13 kBP, the time when the viscous heating peaked (see Figs 2–5 below). It can be noted that viscous dissipation mainly occurs around past ice margins in Laurentia, the North American Cordillera, Fennoscandia and the Barents Sea area where the shear stresses are largest. Due to the constructive interference between the thick Laurentide ice sheets and the thinner Cordilleran ice sheet at 13 kBP, the peak amplitude occurs between their ice margins. However, in the ocean area between the northeast coasts of Canada and Greenland, the destructive interference between ice and water loading results in very low viscous heating. Similarly, constructive interference between the thinner ice sheets in Fennoscandia and the Barents Sea at 13 kBP, results in a smaller local peak between them. Fig. 1 also shows that the magnitude of viscous dissipation for M2 and M3 is more than two times higher than that for M1. The similarity between the results of M2 and M3 at this and all other times indicates that non-linear creep dominates in the composite rheology. However, the presence of an LV upper mantle in M4 results in much higher viscous heating than for M1–M3. The maxima in M1, M2, M3 and M4 are 3.95, 10.14, 10.04 and 22.36 times that of chondritic radiogenic heating, respectively. Figure 1. View largeDownload slide Spatial distribution of the normalized viscous heating for M1–M4, at 13 kBP and depth at 337.5 km. In the subplot for M2, the location of vertical cross-sections in Laurentia and Fennoscandia are indicated by white lines with black dots. Figure 1. View largeDownload slide Spatial distribution of the normalized viscous heating for M1–M4, at 13 kBP and depth at 337.5 km. In the subplot for M2, the location of vertical cross-sections in Laurentia and Fennoscandia are indicated by white lines with black dots. Figure 2. View largeDownload slide Cross-sectional view of viscous heating in Laurentia at different times for M1 and M2. Plots for M3 are not shown because its results are very similar to those of M2. Figure 2. View largeDownload slide Cross-sectional view of viscous heating in Laurentia at different times for M1 and M2. Plots for M3 are not shown because its results are very similar to those of M2. Figure 3. View largeDownload slide Same as in Fig. 2, except for cross-section in Fennoscandia. Figure 3. View largeDownload slide Same as in Fig. 2, except for cross-section in Fennoscandia. Figure 4. View largeDownload slide Cross-sectional view of viscous heating in Laurentia (left-hand panel) and Fennoscandia (right-hand panel) at different times for M4. Figure 4. View largeDownload slide Cross-sectional view of viscous heating in Laurentia (left-hand panel) and Fennoscandia (right-hand panel) at different times for M4. Figure 5. View largeDownload slide (a) The local ice history of the site with maximum viscous heating from 26 thousand years before present time (kBP) to the present. M2 and M3 share the same site which is different from that of M1; (b) the viscous heating dependence on time for that site and (c) the maximum local viscous heating at different times in M1–M3 (see the main text). Figure 5. View largeDownload slide (a) The local ice history of the site with maximum viscous heating from 26 thousand years before present time (kBP) to the present. M2 and M3 share the same site which is different from that of M1; (b) the viscous heating dependence on time for that site and (c) the maximum local viscous heating at different times in M1–M3 (see the main text). To visualize how viscous heating varies with depth, vertical cross-sections in Laurentia and Fennoscandia are shown in Figs 2–4. The locations of these vertical cross-sections are indicated by the two dashed lines in the subplot for M2 in Fig. 1. Figs 2–4 also show how viscous heating varies with time from 15 kBP to 13 and 10 kBP. The spatial-temporal variation in Laurentia for M1 and M2 (similarly for M3) are shown in Fig. 2, while those in Fennoscandia are shown in Fig. 3. Fig. 4 shows the variation in both Laurentia and Fennoscandia for M4. In general, Figs 2–4 show that the location and magnitude of high viscous heat changes with time and in space and the magnitude of viscous heating peaks at 13 kBP. As we shall see in Fig. 5, this is mainly due to the changes of ice thickness and location of the ice centre in time. From Figs 2–4, we see that the peak of viscous heating in M1 generally occurs between 400 and 670 km depth, but in M2 (similarly for M3), the peak under Laurentia is at a shallower depth (Fig. 2). Under Fennoscandia, the locations of the peaks in M1 and M2 are at comparable depth until 13 kBP, after that the peak in M2 is again at a shallower depth than in M1 (Fig. 3). In M4, the presence of the LV upper mantle above a higher viscosity lower mantle pushes the peak to a shallower depth than M1 (Fig. 4). All these are due to the spatial distribution of the von Mises (equivalent) stress τ which is affected by the viscosity structure of the Earth model. Thus, we see that the magnitude of maximum local viscous heating depends on mantle rheology. Fig. 5(a) shows the local ice history at the Earth's surface above the sites with maximum viscous heating for the Earth model under consideration. The time evolution of viscous heating at the sites with maximum heating is shown in Fig. 5(b). At other locations in the Earth, the viscous heating may have a smaller magnitude than the maximum value, but, the value may peak at different times. When all the sites in the mantle are considered, the maximum viscous heating for the whole mantle is shown in Fig. 5(c). Inspection of Figs 5(a) and (b) shows that maximum viscous dissipation occurs near the end of deglaciation at 13 kBP when the rate of melting is fastest. This confirms the finding of the simple ice deglaciation models of Hanyk et al. (2005). In addition, the sharp decline of viscous heating after 13 kBP means that the effect is localized in time, thus the effects of previous glacial cycles are small and can be neglected in our computation. 3.2 Effect of A* and η on the magnitude of viscous heating in uniform mantles In the last subsection, the spatial-temporal variation of viscous heating is shown for a fixed value of A* and η. In this section, we study how the values of A* and η affect the peak value of viscous heating at all times in uniform mantles like M1–M3. We study the range of values for A* from 1 × 10−36 to 1 × 10−33 Pa −3 s−1 and, for η, from 3 × 1020 to 3 × 1022 Pa s as these represent typical values found in uniform mantles. The results are summarized in Table 2. For linear rheology (top row below the labels), it can be seen that the lower the mantle viscosity, the larger is the viscous heating (see discussion of eqs 4a or 4b above). For non-linear rheology (the column on the far right), the largest viscous heating is achieved around A* = 1.11 × 10−34 Pa−3 s−1. This is due to the trade-off between the parameter A* (or A) and τ: for small values of A* (below 1.11 × 10−34 Pa−3 s−1), the relaxation time is long and the magnitude of von Mises stress τ does not change rapidly in time. According to eq. (4a), viscous heating increases with larger value of A* when rheology is non-linear because the effect of stress relaxation is small. However, when A* is above 1.11 × 10−34 Pa−3 s−1, the relaxation time is fast and thus τ decreases rapidly, causing viscous heating to decrease with further increasing of A*. For composite rheology, the trade-off between A* and τ also exist, but the relaxation time of τ is further modified by the value of η. Thus viscous heating no longer peaks at A* around 1.11 × 10−34 Pa−3 s−1 (see Table 2), but at smaller value of A* when the value of η decreases. Table 2 shows that the largest viscous heating within the studied range of parameters is at A* = 1.11 × 10−36 Pa−3 s−1 and when η has the lowest value of 3 × 1020 Pa s (see eq. 4b). 3.3 Effect of mantle compressibility Next, we investigate the effect of material compressibility for linear, non-linear and composite rheology by changing the Poisson's ratio from around 0.5 (incompressible material) to 0.2877 (compressible). The results for models with uniform mantle are summarized in Table 3 which show that compressibility almost doubles the peak viscous heating when the rheology is linear. This is in line with the finding of Hanyk et al. (2005). The reason is that the von Mises stress τ in the upper mantle increases with compressibility. However, from the perspective of strain rate, the reason is less clear: intuitively, with compressibility, volumetric strain rate $${\dot{\varepsilon }_{rr}}^t$$ is greater than zero. However, it is well known (e.g. Chandrasekhar 1981, p. 14, eq. 29) that $$\phi = 2\eta ({\boldsymbol{\dot{\varepsilon }}}_{ij}^2 - \frac{1}{3}\dot{\varepsilon }_{rr}^2)$$ so that ϕ might decrease unless $$\dot{\varepsilon }_{ij}^2$$ increases faster than $$\frac{1}{3}\dot{\varepsilon }_{rr}^2$$, which is apparently the case. Table 3 shows that the effect of compressibility is not as large when the rheology is non-linear or composite because the increase in von Mises stress is smaller. 3.4 Effect of a low viscosity zone From model M4, we saw that an LV upper mantle is able to give a peak value of viscous heating of 22.36 times that of chondritic radiogenic heating. Here, we want to explore the effect of even lower viscosity in the LV zone below the lithosphere for all three types of rheology. The models are modified from VM5a with the addition of an LV zone. The thickness and rheology (A* and η) of the LV zone are listed in Table 4. The linear viscosity in both linear and composite rheology is 3 × 1018 Pa s, that is, two orders of magnitude smaller than that in the upper mantle. The parameter A* is set to 1.11 × 10−34 Pa−3 s−1 in non-linear rheology because this value gives the largest of the peak viscous heating for uniform non-linear rheology models as shown in Table 1. For the same reason, the parameters for composite rheology are taken to be A* = 1.11 × 10−36 Pa−3 s−1 and η = 3 × 1018 Pa s. Results in Table 4 show that the existence of the LV zone can increase the peak value of viscous heating to be over 100 times that of the chondritic radiogenic heating in linear rheology and composite rheology, compared with about 20–30 times in non-linear rheology. This result is also consistent with the predictions of eq. (5). In the following section, we shall investigate whether viscous heating produces significant changes in temperature or surface heat flow. 4 UPPER BOUND ON HEAT FLOW AND TEMPERATURE CHANGE The spatial-temporal evolution of viscous dissipation inside the Earth can be treated as an extra heat source which we denoted as $$H( {\vec{r},t} )$$ and has units of energy generated per unit volume and per unit time. This heat source can raise the local temperature and result in net heat flow according to the heat flow equation (e.g. Turcotte & Schubert 2014),   $$H\left( {\vec{r},t} \right) = \vec{\nabla } \cdot \vec{q}\left( {\vec{r},t} \right) + \rho C\frac{\partial }{{\partial t}}T\left( {\vec{r},t} \right)$$ (5) Here, $$\frac{\partial }{{\partial t}}T( {\vec{r},t} )$$ denotes the rate of temperature change of a volume element with density ρ and specific heat C; $$\vec{q}( {\vec{r},t} )$$ is the heat flux and $$\vec{\nabla } \cdot \vec{q}( {\vec{r},t} )$$ denotes the net heat flow out of the volume element. To complete the computation, note that the heat flux is related to the temperature gradient by Fourier's law of conduction:   $$\vec{q}\left( {\vec{r},t} \right) = - \kappa \vec{\nabla }T\left( {\vec{r},t} \right)$$ (6)where κ is the coefficient of thermal conductivity. In this paper, we are only interested in estimating (i) the upper bound of the induced heat flux (or heat flow) and (ii) the upper bound of the temperature increase due to viscous heating. To estimate the upper bound of the change in heat flux, we take $$\frac{\partial }{{\partial t}}T ( {\vec{r},t} ) = 0$$ in eq. (5) and assume that the heat flow is outwards in the radial direction only. After integrating with respect to the radius, one obtains   $$\delta q \left( {r,\theta ,\varphi ,t} \right) = \mathop \smallint \nolimits_{{\rm{CMB}}}^r H\left( {r',\theta ,\varphi ,t} \right){r^{'2}}dr'/{r^2}$$ (7)where CMB is the core–mantle boundary, and θ and φ are latitude and longitude respectively. To estimate the upper bound of the temperature change, we assume that there is no heat flux in eq. (5) and do the time integration from 26 kBP to the time under consideration (since the heat generated before 26 kBP is very small) to get   $$\delta T \left( {r, \theta, \varphi ,t} \right) = \mathop \smallint \nolimits_{26{\rm{kBP}}}^t H ( {r, theta ,\varphi ,t'} )/( {\rho C} )dt'$$ (8)where C is taken to be 1 kJ kg − 1 °C and uniform in the mantle (e.g. Turcotte & Schubert 2014). Fig. 6 shows the effect of viscous heating on the upper bound of the change in transient heat flux at the surface of the Earth for models M1–M4. The yellow dotted line on M4 indicates the location of the ice margin. Fig. 6 shows that the upper bound of the heat flux at 13 kBP peaks between the edge of the Laurentide and the Cordilleran ice sheets just like in Fig. 1 and the peak magnitude is of the order of 10–30 mW m−2. To put these numbers in perspective, the mean continental heat flux is around 65 mW m−2 (Jaupart & Mareschal 2007) and the mantle heat flux in Canada is about 15 mW m−2 (Jaupart & Mareschal 2007). Thus, the transient heat flux from viscous dissipation due to GIA is not negligible at 13 kBP. An interesting question is whether this extra heat flux near the boundary between the Laurentide and the Cordilleran ice sheets can affect ice stability there. However, since the area of high heat flux is small and the period of high heating is relatively short, the effects on the upper bound of the total heat flow or temperature change are small. This will be demonstrated below. Figure 6. View largeDownload slide Upper bound of perturbed surface heat flux (in mW m−2) at 13 kBP for M1–M4. The yellow dotted line on M4 indicates the location of the ice margin. Figure 6. View largeDownload slide Upper bound of perturbed surface heat flux (in mW m−2) at 13 kBP for M1–M4. The yellow dotted line on M4 indicates the location of the ice margin. How does viscous heating affect the total heat flow out of the Earth? Fig. 7 shows the time evolution of the total heat flow for models M1–M4. This is obtained by integrating the surface heat flux over the surface area of the whole Earth. For M2 and M3, the maximum total heat flow due to viscous heating is at 14 kBP, while for models with linear rheology (i.e. M1 and M4), it is around 10 kBP. The reason why they are at different times is due to the difference in the ηeff and the growth and decay of von Mises stress τ in different models. For M2 and M3, the von Mises stress in the Earth is larger at 14 kBP, while for M1 and M4 it is larger around 10 kBP. For models M1–M4, the total heat flow attained a peak value of around 0.15–0.20 TW, but today the values are less than 0.02 TW. For the models with an LV zone, the peak value of total heat flow is similar to that in M4 although the maximum viscous heating is about five times larger. Since the total heat loss of the Earth today is around 44 TW (Stein 1995), 0.2 TW is less than 0.5 per cent of the current heat loss. Figure 7. View largeDownload slide The time evolution of the upper bound of total perturbed heat flow at the Earth's surface for models M1–M4. Figure 7. View largeDownload slide The time evolution of the upper bound of total perturbed heat flow at the Earth's surface for models M1–M4. Finally, we wish to estimate the upper bound of the temperature change induced by viscous heating for model M4, which has the largest viscous heating in models without LV zone. Fig. 8 shows the upper bound of the temperature change at different depths accumulated up to 13 and 0 kBP (see eq. 8). We can see that the upper bound of the temperature change due to viscous heating is only of the order of 2 × 10−3 °C. However, we can see in Table 4 that the maximum viscous heating for linear and composite rheology models with an LV zone is about five times larger than that for M4, which makes the temperature change almost one order of magnitude larger, that is, 1 × 10−2 °C. Even if there are 30 glacial cycles, then the largest temperature change can only be about 0.3 °C. Thus, the upper bound of the temperature change due to viscous heating is still too small to affect the mantle temperature field, to induce any widespread melting or volcanism or to significantly affect mantle rock properties (including seismic velocities today or mantle viscosity which can potentially affect mantle flow) or affect ice sheet stability. Figure 8. View largeDownload slide The spatial distribution of the upper bound of temperature anomalies (in 10−4 °C) for M4 at depths of 325 and 602.5 km and at 13 and 0 kBP. Figure 8. View largeDownload slide The spatial distribution of the upper bound of temperature anomalies (in 10−4 °C) for M4 at depths of 325 and 602.5 km and at 13 and 0 kBP. 5 CONCLUSIONS We modeled viscous heating in linear, non-linear and composite rheologies using realistic oceans and realistic ice history model ICE-6G. We found that viscous heating is determined by both mantle rheology and ice history. In particular, the magnitude peaks near the end of deglaciation and near the edge of ice sheets. The spatial distribution of viscous heating is also affected by the presence of ocean loading. In general, a small value of η gives the largest viscous heating and for a fixed value of η, larger heating is found in composite rheology than in linear rheology. The viscous heating induced is of the order of 10−9–10−7 W m−3 for a range of reasonable rheology parameters. The upper bound of the transient heat flux due to viscous heating can be as high as about 30 mW m−2 which is not negligible, and is consistent with the finding of Hanyk et al. (2005). However, since the area with heat production is small and the time period of heat production is short, the upper bound of the temperature change due to viscous heating over 30 glacial cycles is of the order of 0.3 °C only. Contrary to the suggestion of Hanyk et al. (2005), we show that even when taking into account non-Newtonian or composite rheology, the upper bound of temperature change due to viscous heating is not large enough to trigger volcanism, affect ice stability or mantle material properties such as seismic velocity and viscosity. ACKNOWLEDGEMENTS We are grateful to Jerry Mitrovica and an anonymous reviewer for their useful comments and suggestions. This work is supported by General Research Fund project17305314 from the Hong Kong Research Grants Council to Patrick Wu. The FE calculation was performed with the ABAQUS package from Hibbitt, Karlsson and Sorensen Inc. This research is conducted in part using the research computing facilities and/or advisory services offered by Information Technology Services, the University of Hong Kong. REFERENCES Chandrasekhar S., 1981. Hydrodynamic and Hydromagnetic Stability , pp. 704, Dover. Gomez N., Pollard D., Holland D., 2015. Sea-level feedback lowers projections of future Antarctic Ice-Sheet mass loss, Nat. Commun. , 6, 8798, doi:10.1038/ncomms9798. https://doi.org/10.1038/ncomms9798 Google Scholar CrossRef Search ADS PubMed  Huybers P., Langmuir C., 2009. Feedback between deglaciation, volcanism, and atmospheric CO2, Earth planet. Sci. Lett. , 286( 3–4), 479– 491. https://doi.org/10.1016/j.epsl.2009.07.014 Google Scholar CrossRef Search ADS   Hanyk L., Matyska C., Yuen D.A., 2005. Short time-scale heating of the Earth's mantle by ice-sheet dynamics, Earth Planet Space , 57( 9), 895– 902. https://doi.org/10.1186/BF03351867 Google Scholar CrossRef Search ADS   Jaupart C., Mareschal J.C., 2007. Heat flow and thermal structure of the Lithosphere-6.05, in Treatise on Geophysics: Crustal and lithosphere dynamics , Vol. 6, pp. 217–253, ed. Schubert G., Elsevier. Jull J., McKenzie D., 1996. The effect of deglaciation on mantle melting beneath Iceland, J. geophys. Res. , 101 ( B10), 21 815– 21 828. https://doi.org/10.1029/96JB01308 Google Scholar CrossRef Search ADS   Karato S.I., Wu P., 1993. Rheology of the upper mantle: a synthesis, Science , 260( 5109), 771– 778. https://doi.org/10.1126/science.260.5109.771 Google Scholar CrossRef Search ADS PubMed  Landau L.D., Lifshitz E.M., 1966. Fluid Mechanics: Landau and Lifshitz: Course of Theoretical Physics , Vol. 6, pp. 50– 51, Elsevier. Nakada M., Yokose H., 1992. Ice age as a trigger of active Quaternary volcanism and tectonism, Tectonophysics , 212( 3–4), 321– 329. https://doi.org/10.1016/0040-1951(92)90298-K Google Scholar CrossRef Search ADS   Pattyn F., 2010, Antarctic subglacial conditions inferred from a hybrid ice sheet/ice stream model, Earth planet. Sci. Lett. , 295( 3–4), 451– 461. https://doi.org/10.1016/j.epsl.2010.04.025 Google Scholar CrossRef Search ADS   Peltier W.R., Argus D.F., Drummond R., 2015. Space geodesy constrains ice age terminal deglaciation: the global ICE-6G_C (VM5a) model, J. geophys. Res. , 120( 1), 450– 487. https://doi.org/10.1002/2014JB011176 Google Scholar CrossRef Search ADS   Ross M.N., Schubert G., 1987. Tidal heating in an internal ocean model of Europa, Nature , 325( 6100), 133– 134. https://doi.org/10.1038/325133a0 Google Scholar CrossRef Search ADS   Ranalli G., 1995. Rheology of the Earth , pp. 413, Springer Science. Sigvaldason G.E., Annertz K., Nilsson M., 1992. Effect of glacier loading/deloading on volcanism: postglacial volcanic production rate of the Dyngjufjöll area, central Iceland, Bull. Volcanol. , 54( 5), 385– 392. https://doi.org/10.1007/BF00312320 Google Scholar CrossRef Search ADS   Segatz M., Spohn T., Ross M.N., Schubert G., 1988. Tidal dissipation, surface heat flow, and figure of viscoelastic models of Io, Icarus , 75( 2), 187– 206. https://doi.org/10.1016/0019-1035(88)90001-2 Google Scholar CrossRef Search ADS   Slater L., Jull M., McKenzie D., Gronvöld K., 1998. Deglaciation effects on mantle melting under Iceland: results from the northern volcanic zone, Earth planet. Sci. Lett. , 164( 1–2), 151– 164. https://doi.org/10.1016/S0012-821X(98)00200-3 Google Scholar CrossRef Search ADS   Stein C.A., 1995. Heat flow of the Earth, in Global Earth Physics: A Handbook of Physical Constants , pp. 144– 158, ed. Ahrens T.J., AGU. Turcotte D.L., Schubert G., 2014. Geodynamics . Cambridge University Press. Uenzelmann-Neben G., Schmidt D.N., Niessen F., Stein R., 2012. Intraplate volcanism off South Greenland: caused by glacial rebound? Geophys. J. Int. , 190( 1), 1– 7. https://doi.org/10.1111/j.1365-246X.2012.05468.x Google Scholar CrossRef Search ADS   van der Wal W., Wu P., Wang H., Sideris M.G., 2010. Sea levels and uplift rate from composite rheology in glacial isostatic adjustment modeling, J. Geodyn. , 50( 1), 38– 48. https://doi.org/10.1016/j.jog.2010.01.006 Google Scholar CrossRef Search ADS   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( 2), 401– 408. https://doi.org/10.1111/j.1365-246X.2004.02338.x Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Effects of mantle rheologies on viscous heating induced by Glacial Isostatic Adjustment

, Volume 213 (1) – Apr 1, 2018
12 pages

/lp/ou_press/effects-of-mantle-rheologies-on-viscous-heating-induced-by-glacial-HKSbjt1w9K
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx535
Publisher site
See Article on Publisher Site

### Abstract

Summary It has been argued that viscous dissipation from mantle flow in response to surface loading during glacial cycles can result in short-term heating and thus trigger transient volcanism or changes in mantle properties, which may in turn affect mantle dynamics. Furthermore, heating near the Earth's surface can also affect the stability of ice sheets. We have studied the magnitude and spatial-temporal distribution of viscous heating induced in the mantle by the realistic ice model ICE-6G and gravitationally consistent ocean loads. Three types of mantle rheologies, including linear, non-linear and composite rheologies are considered to see if non-linear creep can induce larger viscous heating than linear rheology. We used the Coupled-Laplace-Finite-Element model of Glacial Isostatic Adjustment (GIA) to compute the strain, stress and shear heating during a glacial cycle. We also investigated the upper bound of temperature change and surface heat flux change due to viscous heating. We found that maximum viscous heating occurs near the end of deglaciation near the edge of the ice sheet with amplitude as high as 120 times larger than that of the chondritic radioactive heating. The maximum heat flux due to viscous heating can reach 30 mW m−2, but the area with large heat flux is small and the timescale of heating is short. As a result, the upper bound of temperature change due to viscous heating is small. Even if 30 glacial cycles are included, the largest temperature change can be of the order of 0.3 °C. Thus, viscous heating induced by GIA cannot induce volcanism and cannot significantly affect mantle material properties, mantle dynamics nor ice-sheet stability. Loading of the Earth, Europe, North America, Heat generation and transport, Rheology: mantle 1 INTRODUCTION In a fluid, the presence of viscosity, which represents the frictional interaction between neighbouring fluid elements, transforms some kinetic energy of the fluid into heat energy for dissipation (e.g. Landau & Lifshitz 1966). This so-called ‘viscous heating’ can potentially change the temperature and therefore material properties of mantle rocks. If viscous heating is large, then it can be very important because it can induce melting and therefore volcanism. It can also alter seismic properties and thus our interpretation of internal structure of the Earth, or alter mantle viscosity, which may in turn affect mantle dynamics. For example, viscous heating generated by tidal deformation induced by Jupiter may have caused volcanism on the Jovian moon Io (Ross & Schubert 1987; Segatz et al.1988). There is speculation that glacial cycles may have triggered volcanism during the Quaternary (e.g. Nakada & Yokose 1992; Sigvaldason et al.1992; Huybers & Langmuir 2009, Uenzelmann-Neben et al.2012). The link between them is generally attributed to decompression (Jull & McKenzie 1996; Slater et al.1998), but volcanism can also be triggered by viscous heating induced by Glacial Isostatic Adjustment (GIA) of the Earth in response to cycles of ice and ocean loading. Hanyk et al. (2005) modeled viscous heating due to GIA using simple ice load histories on a spherical Earth with linear rheology. They found that viscous heating can generate observable transient heat flux on the surface of the Earth. They also suggested that some degree of volcanism may be associated with viscous heating if the deglaciation period is short, and speculated that higher values of viscous heating can be generated if mantle rheology is non-linear. Indeed, if viscous heating is large enough at the base of the ice sheet, then glacial loading can also affect the stability of ice sheets by basal heating (e.g. Pattyn 2010). On the other hand, if viscous heating leads to lower mantle viscosity near the surface, then this would lead to faster land uplift and that may help to stabilize marine ice sheets (Gomez et al.2015). In this paper, we improve on Hanyk's work by using a realistic ice history model ICE-6G (Peltier et al.2015) together with realistic, self-gravitating oceans. As we shall see, the interaction among ice sheets and the loading of the ocean floor by meltwater affects the spatial distribution of viscous heating. A realistic loading history is also important, as it will eliminate the unrealistic results due to unreasonably fast or slow deglaciation histories. In addition, we consider three types of mantle rheologies in our Earth models: linear, non-linear and composite rheologies to see how large non-linear creep can affect the magnitude of viscous dissipation. Besides studying the effect of viscous heating on surface heat flux, we also study the following questions: what temperature changes inside the Earth are produced by viscous heating induced by glacial cycles? Do these temperature changes significantly affect present-day seismic velocities and thus the interpretation of seismic tomography? Can they affect viscosity and thus mantle flow induced during GIA? In the following sections, we will start with a brief review of the three different types of mantle rheologies, and provide a brief discussion of our model. Then, we will present the results, focusing on the magnitude and spatial and temporal distribution of viscous heating for some simple and other more sophisticated Earth models. After that, we will estimate the effects of viscous heating on the upper bound of heat flux change and see if the upper bound of temperature increase can trigger volcanism and affect ice stability or material properties inside the Earth. 2 THE MODEL 2.1 Mantle rheologies and rate of viscous dissipation Creep experiments on mantle rocks show that both linear and non-linear creep laws operate in the mantle (e.g. Karato & Wu 1993). If both creep laws operate simultaneously, then we have composite rheology. For this case, the creep mechanism with the highest creep rate becomes the dominant creep mechanism. The constitutive relation for these rheologies relates the deviatoric stress σD and the deviatoric strain rate $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$, which are defined as   $${{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}} = {\boldsymbol{\sigma }} - \frac{1}{3}{\sigma _{ii}}{{\bf I}}$$ (1)  $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}} = {\boldsymbol{\dot{\varepsilon }}} - \frac{1}{3}{\dot{\varepsilon }_{ii}}{{\bf I}}$$ (2)where σ and $${\boldsymbol{\dot{\varepsilon }}}$$ are the stress tensor and the strain rate tensor, respectively, I is the identity tensor, σii is the sum of diagonal components of σ and $${\dot{\varepsilon }_{ii}}$$ is the sum of diagonal components of $${\boldsymbol{\dot{\varepsilon }}}$$. In the following, we refer to σDij and $${\dot{\varepsilon }_{Dij}}$$ as the components of σD and $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$, respectively. The relation between σD and $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ for composite, linear and non-linear rheologies can be expressed as (e.g. van der Wal et al.2010):   $${{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}} = \frac{1}{{\frac{1}{{2\eta }} + A{\sigma _E}^{n - 1}}}{{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ (3a) Here, η is the dynamic viscosity for linear creep, A is the non-linear creep parameter determined from shear experiments, and its value depends on temperature, pressure and material properties. σE is the effective deviatoric stress defined by $${\sigma _E} = \sqrt {\frac{1}{2}{\sigma _{Dij}}{\sigma _{Dij}}}$$ and n is the stress exponent with experimental value between 2 and 6. Here, we will take n = 3, as that is a typical value for mantle rocks (Karato & Wu 1993, Ranalli 1995). If A = 0, then eq. (3a) reduces to the relation for linear rheology. If A ≠ 0 and η = ∞, then it becomes the relation for non-linear rheology. It is useful to define the effective viscosity $${\eta _{{\rm{eff}}}} = \frac{1}{2}{( {\frac{1}{{2\eta }} + A{\sigma _E}^{n - 1}} )^{ - 1}}$$, so that eq. (3a) becomes   $${{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}} = 2{\eta _{{\rm{eff}}}}{{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ (3b) For non-linear or composite rheology, the effective viscosity changes with stress level and thus with space and time. For non-linear rheology or composite rheology with a fixed η, a decrease in σE will result in an increase in effective viscosity. The viscous dissipation rate ϕ (hereafter called ‘viscous heating’) is given by (e.g. eq. 7 in Hanyk et al.2005):   $$\phi = {{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}}{:}\, {{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}} = \frac{1}{{2{{\boldsymbol{\eta }}_{{{\bf eff}}}}}}{{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}}{:}\,{{\boldsymbol{\sigma\! }}_{\boldsymbol{D}}}$$ (4a) This shows that shear heating decreases if the effective viscosity increases provided that σD does not change rapidly. For surface loading problems like GIA, σD is determined by the changing surface loads, although σD is also affected by stress relaxation where the decay time constant is proportional to the effective viscosity. In terms of von Mises stress $$\tau = \sqrt 3 {\sigma _E}$$, the creep parameter A and viscosity η, eq. (4a) can be written as:   $$\phi = \frac{2}{3} \left( {\frac{1}{{2\eta }} + \frac{1}{3}A{\tau ^2}} \right){\tau ^2}$$ (4b)which clearly shows the strong dependence of ϕ on von Mises stress τ. Eq. (4a) can also be written as: $$\phi = 2{{\boldsymbol{\eta }}_{{{\bf eff}}}}{{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}{:}\,{{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ which shows that viscous heating would increase with larger effective viscosity provided that the strain rate $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ is constant. For our problem, $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ is not constant, but $${{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}{:}\,{{\boldsymbol{\dot{\varepsilon\! }}}_{\boldsymbol{D}}}$$ is determined by (σD/ηeff)2. Thus, eqs (4a) or (4b) will be used for our discussion below. 2.2 The GIA model The Coupled Laplace-Finite-Element (CLFE) GIA model of Wu (2004), later modified by van der Wal et al. (2010), is used in this study. Finite-element grids with different (0.5° and 2°) spatial resolution have been used for the computations. It is found that 2° resolution is adequate for our purpose. The inputs of the GIA model are elastic and creep parameters for the Earth model, and ice loading history. Unlike the simple ice sheet model used in Hanyk et al. (2005), the ice model used here is the realistic global ICE6G model of Peltier et al. (2015). Since ICE6G provides ice thickness history from 26 kBP to the present, we assume that the ice thickness increased linearly from zero at 108 kBP to the ice thickness at 26 kBP. Also, the ice was taken from the water in the oceans and meltwater returned to the oceans. The self-consistent sea level equation is solved for realistic oceans (Wu 2004). Since the effects of time-dependent coastline and rotational feedback are small on the von Mises stress τ and shear heating, they have not been included. As we shall see below, the localization of viscous heating in time makes the consideration of previous glacial cycles unnecessary. For the Earth models, the elastic parameters for models M1–M3 are the same as those in van der Wal et al. (2010), and for M4, they are the same as that in model VM5a of Peltier et al. (2015). The rheological parameters are described in Table 1. Here, parameter A* is the creep parameter determined from uniaxial experiments and is related to A, the creep parameter from shear experiments by $${A^*} = 2A/\sqrt {{3^{n + 1}}}$$ (van der Wal et al.2010). In Table 1, M1 and M4 are linear rheological models. M1 has a uniform viscosity of 3 × 1021 Pa s, while in M4, the viscosities in the upper mantle (100–670 km depth), shallow lower mantle (670–1271 km depth) and deep lower mantle (1271–2891 km depth) are 5 × 1020, 1.6 × 1021 and 3 × 1021 Pa s respectively. M4 is considered because model VM5a is used for the construction of ICE6G (Peltier et al.2015), and so these should be used together. The rheological parameters in models M2 and M3 are uniform in the mantle. M2 has a non-linear rheology, while model M3 has a composite rheology. In other models, the effects of compressibility and a low viscosity (LV) zone below the lithosphere are also studied (see Tables 2–4) as they can also affect the magnitude of viscous heating. Table 1. Rheological parameters of four Earth models.   M1  M2  M3  M4  A*(Pa − 3 · s − 1)  0  1.11 × 10−34  1.11 × 10−34  0  η (Pa s)  3 × 1021  0  3 × 1021  VM5a  n  1  3  3  1    M1  M2  M3  M4  A*(Pa − 3 · s − 1)  0  1.11 × 10−34  1.11 × 10−34  0  η (Pa s)  3 × 1021  0  3 × 1021  VM5a  n  1  3  3  1  View Large Table 2. The maximum local viscous heating of all time for uniform mantle models with various rheology (A* and η). ϕMax  η(Pa · s)      3.00 × 1020  3.00 × 1021  3.00 × 1022  Non-linear  A* (Pa − 3 · s − 1)  Linear  11.64  3.95  0.54      1.11 × 10−36  11.64  6.14  2.73  2.23    1.11 × 10−35  11.45  11.45  10.24  9.99    1.11 × 10−34  9.58  10.04  10.12  10.14    1.11 × 10−33  6.73  6.55  6.54  6.53  ϕMax  η(Pa · s)      3.00 × 1020  3.00 × 1021  3.00 × 1022  Non-linear  A* (Pa − 3 · s − 1)  Linear  11.64  3.95  0.54      1.11 × 10−36  11.64  6.14  2.73  2.23    1.11 × 10−35  11.45  11.45  10.24  9.99    1.11 × 10−34  9.58  10.04  10.12  10.14    1.11 × 10−33  6.73  6.55  6.54  6.53  View Large Table 3. The maximum local viscous heating of all epochs for two values of Poisson's ratio (i.e. compressibility) with various rheology (A* [Pa−3 s−1] and η [Pa s]). ϕMax  A* = 0 η = 3 × 1021  A* = 1.11 × 10−34 (non − linear)  A* = 1.11 × 10−34η = 3 × 1021  Poisson's ratio  0.4900  3.95  10.14  10.04    0.2877  6.18  10.59  10.40  ϕMax  A* = 0 η = 3 × 1021  A* = 1.11 × 10−34 (non − linear)  A* = 1.11 × 10−34η = 3 × 1021  Poisson's ratio  0.4900  3.95  10.14  10.04    0.2877  6.18  10.59  10.40  View Large Table 4. The maximum local viscous heat of all time in models with an LV zone below the lithosphere. Units: A∗ in [Pa−3 s−1] and η in [Pa s]. ϕMax  A*=0 η = 3 × 1018  A* = 1.11 × 10−34 (non − linear)  A* = 1.11 × 10−36η = 3 × 1018  LV zone thickness (km)  40  102.18  26.01  102.19    100  118.80  30.52  118.81  ϕMax  A*=0 η = 3 × 1018  A* = 1.11 × 10−34 (non − linear)  A* = 1.11 × 10−36η = 3 × 1018  LV zone thickness (km)  40  102.18  26.01  102.19    100  118.80  30.52  118.81  View Large The outputs of the CLFE model contain spatial temporal evolution of the displacements and state of stress throughout the mantle. From the stress output at any time step, the local viscous heating for each element can be computed using eqs (4a) or (4b). In this paper, viscous heating is normalized by the chondritic radiogenic heating of 3 × 10−9 W m−3 (Hanyk et al.2005). 3 RESULTS 3.1 Spatial-temporal distribution of viscous heating for M1–M4 In Fig. 1, the spatial distribution of normalized viscous heating for models M1–M4 is shown at 13 kBP, the time when the viscous heating peaked (see Figs 2–5 below). It can be noted that viscous dissipation mainly occurs around past ice margins in Laurentia, the North American Cordillera, Fennoscandia and the Barents Sea area where the shear stresses are largest. Due to the constructive interference between the thick Laurentide ice sheets and the thinner Cordilleran ice sheet at 13 kBP, the peak amplitude occurs between their ice margins. However, in the ocean area between the northeast coasts of Canada and Greenland, the destructive interference between ice and water loading results in very low viscous heating. Similarly, constructive interference between the thinner ice sheets in Fennoscandia and the Barents Sea at 13 kBP, results in a smaller local peak between them. Fig. 1 also shows that the magnitude of viscous dissipation for M2 and M3 is more than two times higher than that for M1. The similarity between the results of M2 and M3 at this and all other times indicates that non-linear creep dominates in the composite rheology. However, the presence of an LV upper mantle in M4 results in much higher viscous heating than for M1–M3. The maxima in M1, M2, M3 and M4 are 3.95, 10.14, 10.04 and 22.36 times that of chondritic radiogenic heating, respectively. Figure 1. View largeDownload slide Spatial distribution of the normalized viscous heating for M1–M4, at 13 kBP and depth at 337.5 km. In the subplot for M2, the location of vertical cross-sections in Laurentia and Fennoscandia are indicated by white lines with black dots. Figure 1. View largeDownload slide Spatial distribution of the normalized viscous heating for M1–M4, at 13 kBP and depth at 337.5 km. In the subplot for M2, the location of vertical cross-sections in Laurentia and Fennoscandia are indicated by white lines with black dots. Figure 2. View largeDownload slide Cross-sectional view of viscous heating in Laurentia at different times for M1 and M2. Plots for M3 are not shown because its results are very similar to those of M2. Figure 2. View largeDownload slide Cross-sectional view of viscous heating in Laurentia at different times for M1 and M2. Plots for M3 are not shown because its results are very similar to those of M2. Figure 3. View largeDownload slide Same as in Fig. 2, except for cross-section in Fennoscandia. Figure 3. View largeDownload slide Same as in Fig. 2, except for cross-section in Fennoscandia. Figure 4. View largeDownload slide Cross-sectional view of viscous heating in Laurentia (left-hand panel) and Fennoscandia (right-hand panel) at different times for M4. Figure 4. View largeDownload slide Cross-sectional view of viscous heating in Laurentia (left-hand panel) and Fennoscandia (right-hand panel) at different times for M4. Figure 5. View largeDownload slide (a) The local ice history of the site with maximum viscous heating from 26 thousand years before present time (kBP) to the present. M2 and M3 share the same site which is different from that of M1; (b) the viscous heating dependence on time for that site and (c) the maximum local viscous heating at different times in M1–M3 (see the main text). Figure 5. View largeDownload slide (a) The local ice history of the site with maximum viscous heating from 26 thousand years before present time (kBP) to the present. M2 and M3 share the same site which is different from that of M1; (b) the viscous heating dependence on time for that site and (c) the maximum local viscous heating at different times in M1–M3 (see the main text). To visualize how viscous heating varies with depth, vertical cross-sections in Laurentia and Fennoscandia are shown in Figs 2–4. The locations of these vertical cross-sections are indicated by the two dashed lines in the subplot for M2 in Fig. 1. Figs 2–4 also show how viscous heating varies with time from 15 kBP to 13 and 10 kBP. The spatial-temporal variation in Laurentia for M1 and M2 (similarly for M3) are shown in Fig. 2, while those in Fennoscandia are shown in Fig. 3. Fig. 4 shows the variation in both Laurentia and Fennoscandia for M4. In general, Figs 2–4 show that the location and magnitude of high viscous heat changes with time and in space and the magnitude of viscous heating peaks at 13 kBP. As we shall see in Fig. 5, this is mainly due to the changes of ice thickness and location of the ice centre in time. From Figs 2–4, we see that the peak of viscous heating in M1 generally occurs between 400 and 670 km depth, but in M2 (similarly for M3), the peak under Laurentia is at a shallower depth (Fig. 2). Under Fennoscandia, the locations of the peaks in M1 and M2 are at comparable depth until 13 kBP, after that the peak in M2 is again at a shallower depth than in M1 (Fig. 3). In M4, the presence of the LV upper mantle above a higher viscosity lower mantle pushes the peak to a shallower depth than M1 (Fig. 4). All these are due to the spatial distribution of the von Mises (equivalent) stress τ which is affected by the viscosity structure of the Earth model. Thus, we see that the magnitude of maximum local viscous heating depends on mantle rheology. Fig. 5(a) shows the local ice history at the Earth's surface above the sites with maximum viscous heating for the Earth model under consideration. The time evolution of viscous heating at the sites with maximum heating is shown in Fig. 5(b). At other locations in the Earth, the viscous heating may have a smaller magnitude than the maximum value, but, the value may peak at different times. When all the sites in the mantle are considered, the maximum viscous heating for the whole mantle is shown in Fig. 5(c). Inspection of Figs 5(a) and (b) shows that maximum viscous dissipation occurs near the end of deglaciation at 13 kBP when the rate of melting is fastest. This confirms the finding of the simple ice deglaciation models of Hanyk et al. (2005). In addition, the sharp decline of viscous heating after 13 kBP means that the effect is localized in time, thus the effects of previous glacial cycles are small and can be neglected in our computation. 3.2 Effect of A* and η on the magnitude of viscous heating in uniform mantles In the last subsection, the spatial-temporal variation of viscous heating is shown for a fixed value of A* and η. In this section, we study how the values of A* and η affect the peak value of viscous heating at all times in uniform mantles like M1–M3. We study the range of values for A* from 1 × 10−36 to 1 × 10−33 Pa −3 s−1 and, for η, from 3 × 1020 to 3 × 1022 Pa s as these represent typical values found in uniform mantles. The results are summarized in Table 2. For linear rheology (top row below the labels), it can be seen that the lower the mantle viscosity, the larger is the viscous heating (see discussion of eqs 4a or 4b above). For non-linear rheology (the column on the far right), the largest viscous heating is achieved around A* = 1.11 × 10−34 Pa−3 s−1. This is due to the trade-off between the parameter A* (or A) and τ: for small values of A* (below 1.11 × 10−34 Pa−3 s−1), the relaxation time is long and the magnitude of von Mises stress τ does not change rapidly in time. According to eq. (4a), viscous heating increases with larger value of A* when rheology is non-linear because the effect of stress relaxation is small. However, when A* is above 1.11 × 10−34 Pa−3 s−1, the relaxation time is fast and thus τ decreases rapidly, causing viscous heating to decrease with further increasing of A*. For composite rheology, the trade-off between A* and τ also exist, but the relaxation time of τ is further modified by the value of η. Thus viscous heating no longer peaks at A* around 1.11 × 10−34 Pa−3 s−1 (see Table 2), but at smaller value of A* when the value of η decreases. Table 2 shows that the largest viscous heating within the studied range of parameters is at A* = 1.11 × 10−36 Pa−3 s−1 and when η has the lowest value of 3 × 1020 Pa s (see eq. 4b). 3.3 Effect of mantle compressibility Next, we investigate the effect of material compressibility for linear, non-linear and composite rheology by changing the Poisson's ratio from around 0.5 (incompressible material) to 0.2877 (compressible). The results for models with uniform mantle are summarized in Table 3 which show that compressibility almost doubles the peak viscous heating when the rheology is linear. This is in line with the finding of Hanyk et al. (2005). The reason is that the von Mises stress τ in the upper mantle increases with compressibility. However, from the perspective of strain rate, the reason is less clear: intuitively, with compressibility, volumetric strain rate $${\dot{\varepsilon }_{rr}}^t$$ is greater than zero. However, it is well known (e.g. Chandrasekhar 1981, p. 14, eq. 29) that $$\phi = 2\eta ({\boldsymbol{\dot{\varepsilon }}}_{ij}^2 - \frac{1}{3}\dot{\varepsilon }_{rr}^2)$$ so that ϕ might decrease unless $$\dot{\varepsilon }_{ij}^2$$ increases faster than $$\frac{1}{3}\dot{\varepsilon }_{rr}^2$$, which is apparently the case. Table 3 shows that the effect of compressibility is not as large when the rheology is non-linear or composite because the increase in von Mises stress is smaller. 3.4 Effect of a low viscosity zone From model M4, we saw that an LV upper mantle is able to give a peak value of viscous heating of 22.36 times that of chondritic radiogenic heating. Here, we want to explore the effect of even lower viscosity in the LV zone below the lithosphere for all three types of rheology. The models are modified from VM5a with the addition of an LV zone. The thickness and rheology (A* and η) of the LV zone are listed in Table 4. The linear viscosity in both linear and composite rheology is 3 × 1018 Pa s, that is, two orders of magnitude smaller than that in the upper mantle. The parameter A* is set to 1.11 × 10−34 Pa−3 s−1 in non-linear rheology because this value gives the largest of the peak viscous heating for uniform non-linear rheology models as shown in Table 1. For the same reason, the parameters for composite rheology are taken to be A* = 1.11 × 10−36 Pa−3 s−1 and η = 3 × 1018 Pa s. Results in Table 4 show that the existence of the LV zone can increase the peak value of viscous heating to be over 100 times that of the chondritic radiogenic heating in linear rheology and composite rheology, compared with about 20–30 times in non-linear rheology. This result is also consistent with the predictions of eq. (5). In the following section, we shall investigate whether viscous heating produces significant changes in temperature or surface heat flow. 4 UPPER BOUND ON HEAT FLOW AND TEMPERATURE CHANGE The spatial-temporal evolution of viscous dissipation inside the Earth can be treated as an extra heat source which we denoted as $$H( {\vec{r},t} )$$ and has units of energy generated per unit volume and per unit time. This heat source can raise the local temperature and result in net heat flow according to the heat flow equation (e.g. Turcotte & Schubert 2014),   $$H\left( {\vec{r},t} \right) = \vec{\nabla } \cdot \vec{q}\left( {\vec{r},t} \right) + \rho C\frac{\partial }{{\partial t}}T\left( {\vec{r},t} \right)$$ (5) Here, $$\frac{\partial }{{\partial t}}T( {\vec{r},t} )$$ denotes the rate of temperature change of a volume element with density ρ and specific heat C; $$\vec{q}( {\vec{r},t} )$$ is the heat flux and $$\vec{\nabla } \cdot \vec{q}( {\vec{r},t} )$$ denotes the net heat flow out of the volume element. To complete the computation, note that the heat flux is related to the temperature gradient by Fourier's law of conduction:   $$\vec{q}\left( {\vec{r},t} \right) = - \kappa \vec{\nabla }T\left( {\vec{r},t} \right)$$ (6)where κ is the coefficient of thermal conductivity. In this paper, we are only interested in estimating (i) the upper bound of the induced heat flux (or heat flow) and (ii) the upper bound of the temperature increase due to viscous heating. To estimate the upper bound of the change in heat flux, we take $$\frac{\partial }{{\partial t}}T ( {\vec{r},t} ) = 0$$ in eq. (5) and assume that the heat flow is outwards in the radial direction only. After integrating with respect to the radius, one obtains   $$\delta q \left( {r,\theta ,\varphi ,t} \right) = \mathop \smallint \nolimits_{{\rm{CMB}}}^r H\left( {r',\theta ,\varphi ,t} \right){r^{'2}}dr'/{r^2}$$ (7)where CMB is the core–mantle boundary, and θ and φ are latitude and longitude respectively. To estimate the upper bound of the temperature change, we assume that there is no heat flux in eq. (5) and do the time integration from 26 kBP to the time under consideration (since the heat generated before 26 kBP is very small) to get   $$\delta T \left( {r, \theta, \varphi ,t} \right) = \mathop \smallint \nolimits_{26{\rm{kBP}}}^t H ( {r, theta ,\varphi ,t'} )/( {\rho C} )dt'$$ (8)where C is taken to be 1 kJ kg − 1 °C and uniform in the mantle (e.g. Turcotte & Schubert 2014). Fig. 6 shows the effect of viscous heating on the upper bound of the change in transient heat flux at the surface of the Earth for models M1–M4. The yellow dotted line on M4 indicates the location of the ice margin. Fig. 6 shows that the upper bound of the heat flux at 13 kBP peaks between the edge of the Laurentide and the Cordilleran ice sheets just like in Fig. 1 and the peak magnitude is of the order of 10–30 mW m−2. To put these numbers in perspective, the mean continental heat flux is around 65 mW m−2 (Jaupart & Mareschal 2007) and the mantle heat flux in Canada is about 15 mW m−2 (Jaupart & Mareschal 2007). Thus, the transient heat flux from viscous dissipation due to GIA is not negligible at 13 kBP. An interesting question is whether this extra heat flux near the boundary between the Laurentide and the Cordilleran ice sheets can affect ice stability there. However, since the area of high heat flux is small and the period of high heating is relatively short, the effects on the upper bound of the total heat flow or temperature change are small. This will be demonstrated below. Figure 6. View largeDownload slide Upper bound of perturbed surface heat flux (in mW m−2) at 13 kBP for M1–M4. The yellow dotted line on M4 indicates the location of the ice margin. Figure 6. View largeDownload slide Upper bound of perturbed surface heat flux (in mW m−2) at 13 kBP for M1–M4. The yellow dotted line on M4 indicates the location of the ice margin. How does viscous heating affect the total heat flow out of the Earth? Fig. 7 shows the time evolution of the total heat flow for models M1–M4. This is obtained by integrating the surface heat flux over the surface area of the whole Earth. For M2 and M3, the maximum total heat flow due to viscous heating is at 14 kBP, while for models with linear rheology (i.e. M1 and M4), it is around 10 kBP. The reason why they are at different times is due to the difference in the ηeff and the growth and decay of von Mises stress τ in different models. For M2 and M3, the von Mises stress in the Earth is larger at 14 kBP, while for M1 and M4 it is larger around 10 kBP. For models M1–M4, the total heat flow attained a peak value of around 0.15–0.20 TW, but today the values are less than 0.02 TW. For the models with an LV zone, the peak value of total heat flow is similar to that in M4 although the maximum viscous heating is about five times larger. Since the total heat loss of the Earth today is around 44 TW (Stein 1995), 0.2 TW is less than 0.5 per cent of the current heat loss. Figure 7. View largeDownload slide The time evolution of the upper bound of total perturbed heat flow at the Earth's surface for models M1–M4. Figure 7. View largeDownload slide The time evolution of the upper bound of total perturbed heat flow at the Earth's surface for models M1–M4. Finally, we wish to estimate the upper bound of the temperature change induced by viscous heating for model M4, which has the largest viscous heating in models without LV zone. Fig. 8 shows the upper bound of the temperature change at different depths accumulated up to 13 and 0 kBP (see eq. 8). We can see that the upper bound of the temperature change due to viscous heating is only of the order of 2 × 10−3 °C. However, we can see in Table 4 that the maximum viscous heating for linear and composite rheology models with an LV zone is about five times larger than that for M4, which makes the temperature change almost one order of magnitude larger, that is, 1 × 10−2 °C. Even if there are 30 glacial cycles, then the largest temperature change can only be about 0.3 °C. Thus, the upper bound of the temperature change due to viscous heating is still too small to affect the mantle temperature field, to induce any widespread melting or volcanism or to significantly affect mantle rock properties (including seismic velocities today or mantle viscosity which can potentially affect mantle flow) or affect ice sheet stability. Figure 8. View largeDownload slide The spatial distribution of the upper bound of temperature anomalies (in 10−4 °C) for M4 at depths of 325 and 602.5 km and at 13 and 0 kBP. Figure 8. View largeDownload slide The spatial distribution of the upper bound of temperature anomalies (in 10−4 °C) for M4 at depths of 325 and 602.5 km and at 13 and 0 kBP. 5 CONCLUSIONS We modeled viscous heating in linear, non-linear and composite rheologies using realistic oceans and realistic ice history model ICE-6G. We found that viscous heating is determined by both mantle rheology and ice history. In particular, the magnitude peaks near the end of deglaciation and near the edge of ice sheets. The spatial distribution of viscous heating is also affected by the presence of ocean loading. In general, a small value of η gives the largest viscous heating and for a fixed value of η, larger heating is found in composite rheology than in linear rheology. The viscous heating induced is of the order of 10−9–10−7 W m−3 for a range of reasonable rheology parameters. The upper bound of the transient heat flux due to viscous heating can be as high as about 30 mW m−2 which is not negligible, and is consistent with the finding of Hanyk et al. (2005). However, since the area with heat production is small and the time period of heat production is short, the upper bound of the temperature change due to viscous heating over 30 glacial cycles is of the order of 0.3 °C only. Contrary to the suggestion of Hanyk et al. (2005), we show that even when taking into account non-Newtonian or composite rheology, the upper bound of temperature change due to viscous heating is not large enough to trigger volcanism, affect ice stability or mantle material properties such as seismic velocity and viscosity. ACKNOWLEDGEMENTS We are grateful to Jerry Mitrovica and an anonymous reviewer for their useful comments and suggestions. This work is supported by General Research Fund project17305314 from the Hong Kong Research Grants Council to Patrick Wu. The FE calculation was performed with the ABAQUS package from Hibbitt, Karlsson and Sorensen Inc. This research is conducted in part using the research computing facilities and/or advisory services offered by Information Technology Services, the University of Hong Kong. REFERENCES Chandrasekhar S., 1981. Hydrodynamic and Hydromagnetic Stability , pp. 704, Dover. Gomez N., Pollard D., Holland D., 2015. Sea-level feedback lowers projections of future Antarctic Ice-Sheet mass loss, Nat. Commun. , 6, 8798, doi:10.1038/ncomms9798. https://doi.org/10.1038/ncomms9798 Google Scholar CrossRef Search ADS PubMed  Huybers P., Langmuir C., 2009. Feedback between deglaciation, volcanism, and atmospheric CO2, Earth planet. Sci. Lett. , 286( 3–4), 479– 491. https://doi.org/10.1016/j.epsl.2009.07.014 Google Scholar CrossRef Search ADS   Hanyk L., Matyska C., Yuen D.A., 2005. Short time-scale heating of the Earth's mantle by ice-sheet dynamics, Earth Planet Space , 57( 9), 895– 902. https://doi.org/10.1186/BF03351867 Google Scholar CrossRef Search ADS   Jaupart C., Mareschal J.C., 2007. Heat flow and thermal structure of the Lithosphere-6.05, in Treatise on Geophysics: Crustal and lithosphere dynamics , Vol. 6, pp. 217–253, ed. Schubert G., Elsevier. Jull J., McKenzie D., 1996. The effect of deglaciation on mantle melting beneath Iceland, J. geophys. Res. , 101 ( B10), 21 815– 21 828. https://doi.org/10.1029/96JB01308 Google Scholar CrossRef Search ADS   Karato S.I., Wu P., 1993. Rheology of the upper mantle: a synthesis, Science , 260( 5109), 771– 778. https://doi.org/10.1126/science.260.5109.771 Google Scholar CrossRef Search ADS PubMed  Landau L.D., Lifshitz E.M., 1966. Fluid Mechanics: Landau and Lifshitz: Course of Theoretical Physics , Vol. 6, pp. 50– 51, Elsevier. Nakada M., Yokose H., 1992. Ice age as a trigger of active Quaternary volcanism and tectonism, Tectonophysics , 212( 3–4), 321– 329. https://doi.org/10.1016/0040-1951(92)90298-K Google Scholar CrossRef Search ADS   Pattyn F., 2010, Antarctic subglacial conditions inferred from a hybrid ice sheet/ice stream model, Earth planet. Sci. Lett. , 295( 3–4), 451– 461. https://doi.org/10.1016/j.epsl.2010.04.025 Google Scholar CrossRef Search ADS   Peltier W.R., Argus D.F., Drummond R., 2015. Space geodesy constrains ice age terminal deglaciation: the global ICE-6G_C (VM5a) model, J. geophys. Res. , 120( 1), 450– 487. https://doi.org/10.1002/2014JB011176 Google Scholar CrossRef Search ADS   Ross M.N., Schubert G., 1987. Tidal heating in an internal ocean model of Europa, Nature , 325( 6100), 133– 134. https://doi.org/10.1038/325133a0 Google Scholar CrossRef Search ADS   Ranalli G., 1995. Rheology of the Earth , pp. 413, Springer Science. Sigvaldason G.E., Annertz K., Nilsson M., 1992. Effect of glacier loading/deloading on volcanism: postglacial volcanic production rate of the Dyngjufjöll area, central Iceland, Bull. Volcanol. , 54( 5), 385– 392. https://doi.org/10.1007/BF00312320 Google Scholar CrossRef Search ADS   Segatz M., Spohn T., Ross M.N., Schubert G., 1988. Tidal dissipation, surface heat flow, and figure of viscoelastic models of Io, Icarus , 75( 2), 187– 206. https://doi.org/10.1016/0019-1035(88)90001-2 Google Scholar CrossRef Search ADS   Slater L., Jull M., McKenzie D., Gronvöld K., 1998. Deglaciation effects on mantle melting under Iceland: results from the northern volcanic zone, Earth planet. Sci. Lett. , 164( 1–2), 151– 164. https://doi.org/10.1016/S0012-821X(98)00200-3 Google Scholar CrossRef Search ADS   Stein C.A., 1995. Heat flow of the Earth, in Global Earth Physics: A Handbook of Physical Constants , pp. 144– 158, ed. Ahrens T.J., AGU. Turcotte D.L., Schubert G., 2014. Geodynamics . Cambridge University Press. Uenzelmann-Neben G., Schmidt D.N., Niessen F., Stein R., 2012. Intraplate volcanism off South Greenland: caused by glacial rebound? Geophys. J. Int. , 190( 1), 1– 7. https://doi.org/10.1111/j.1365-246X.2012.05468.x Google Scholar CrossRef Search ADS   van der Wal W., Wu P., Wang H., Sideris M.G., 2010. Sea levels and uplift rate from composite rheology in glacial isostatic adjustment modeling, J. Geodyn. , 50( 1), 38– 48. https://doi.org/10.1016/j.jog.2010.01.006 Google Scholar CrossRef Search ADS   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( 2), 401– 408. https://doi.org/10.1111/j.1365-246X.2004.02338.x Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations