Incorporating Deformation Energetics in Long-Term Tectonic Modeling
aa r X i v : . [ phy s i c s . g e o - ph ] D ec Incorporating Deformation Energetics in Long-Term TectonicModeling
Sabber Ahamed a, ∗ , Eunseo Choi a a Center for Earthquake Research and Information, The University of Memphis, 3890 Central Ave.,Memphis, TN 38152, U.S.A.
Abstract
The deformation-related energy budget is usually considered in the simplest formor even completely omitted from the energy balance equation. We derive a fullenergy balance equation that accounts not only for heat energy but also for me-chanical (elastic, plastic and viscous) work. The derived equation is implementedin DES3D, an unstructured finite element solver for long-term tectonic deforma-tion. We verify the implementation by comparing numerical solutions to the corre-sponding semi-analytic solutions in three benchmarks extended from the classicaloedometer test. Two of the benchmarks are designed to evaluate the temperaturechange in a Mohr-Coulomb elasto-plastic square governed by a simplified equa-tion involving plastic power only and by the full temperature evolution equation,respectively. The third benchmark di ff ers in that it computes thermal stresses as-sociated with a prescribed uniform temperature increase. All the solutions fromDES3D show relative errors less than 0.1 %. We also investigate the long-terme ff ects of deformation energetics on the evolution of large o ff set normal faults. Wefind that the models considering the full energy balance equation tend to producemore secondary faults and an elongated core complex. Our results for the normal ∗ Corresponding author
Email address: [email protected] (Sabber Ahamed)
Preprint submitted to Tectonophysics December 27, 2018 ault system confirm that persistent inelastic deformation has significant impact onthe long-term evolution of faults, motivating further exploration of the role of thefull energy balance equation in other geodynamic systems.
Keywords: plastic power, strain-softening plasticity, thermodynamic principles,thermal stress
1. Introduction
The energy conservation principle can describe a wider range of geologicaland geodynamic phenomena when it takes into account energies involved withthan the usual form involving only heat advection and di ff usion because the gen-eral form can account for the energetics of deformation, which is often impor-tant in complicated phenomena. For instance, Hunt and Wadee (1991) show thatnon-periodic localized folding can be viewed as a superposition of folding modesthat correspond to multiple local minima of the non-convex energy function andthe non-convexity originates from deformational contribution to the system’s en-ergy budget. On a similar note, Hobbs et al. (2011) propose that the feedback be-tween shear heating and temperature-dependent viscosity can explain non-periodicand non-symmetric folding occurring in layers with small viscosity contrast, Theclassical Biot’s theory (e.g., Biot, 1961) predicts that folding would not occur insuch a configuration. Influences of energy dissipation have been considered in thelithospheric scale as well. Regenauer-Lieb et al. (2006) consider the feedback be-tween the energy dissipated due to inelastic deformation and changes in viscousand plastic material properties due to temperature changes resulting from heat con-verted from the dissipated energy. They find that the two-way feedback processcan make the brittle-ductile transition (BDT) zone weaker than other parts of litho-sphere although the classical strength envelopes predict the BDT zone of litho-sphere to be the strongest (Ranalli and Murphy, 1987; Goetze and Evans, 1979;2race and Kohlstedt, 1980). Energy dissipated in the form of shear heating isshown to promote a necking in the subducting slab and ultimately lead to slab de-tachment (Gerya et al., 2004). In the whole-mantle scale, Yuen et al. (1987) pointout that feedback between rheology and dissipative energy in mantle convectionis potentially an important mechanism that can warm up the mantle by severalhundred degrees above the incompressible profile. Ita and King (1994) show thatdissipative heating in the energy balance can significantly facilitate the vertical flowacross the 660-km phase boundary.However, computational long-term tectonic modeling, concerned about thelong-term evolution of geological structures of various scales, has yet to fully em-brace deformation energetics. Problems are centered around the fact that simpli-fications commonly made in this type of modeling preclude consistent thermo-mechanical coupling. For instance, energy balance is considered only in terms ofheat advection and di ff usion while thermo-mechanical feedback is realized onlythrough temperature-dependent viscosity or shear heating.More specifically, we identify three elements in kinematics and constitutivemodels that need to be incorporated into a thermodyanmic framework for long-term tectonic models. Firstly, we note that the elastic or plastic deformations arefrequently assumed to be incompressible (e.g., Regenauer-Lieb and Yuen, 2003;Regenauer-Lieb et al., 2006, 2008; Connolly, 2009; Hobbs et al., 2011) althoughthis assumption is neither required nor well-justified. Volumetric strain can have asignificant e ff ect on energy budget (Hunsche, 1991; Zinoviev and Ermakov, 1994)and is non-negligible during brittle deformations (e.g. Brace et al., 1966) and phasetransformations Hyndman and Peacock (2003); Het´enyi et al. (2011). Secondly,thermal stresses are often ignored even though they can be a significant source oftransient stresses and associated deformation (Choi et al., 2008; Choi and Gurnis,2008; Korenaga, 2007). A simple back-of-the-envelope calculation shows that a3emperature change of 100 K in a perfectly confined rock body with a bulk modulusof 30 GPa can generate up to 90 MPa thermal stresses when the thermal expansioncoe ffi cient is 3 × − K − . Although transient, thermal stresses of such a magni-tude might be su ffi cient for driving permanent changes in state variables such aselastic damage or plastic strain under non-linear rheologies. Thirdly, strain weak-ening plasticity is sometimes considered as inconsistent with thermodynamic prin-ciples (e.g. Regenauer-Lieb et al., 2006). However, frictional materials do show re-duction in overall strength with continued loading (e.g. Read and Hegemier, 1984;Borja et al., 2000), making it necessary to consider the softening behavior in tec-tonic models concerned about brittle behaviors of rocks. In fact, the strain softeningis perfectly legitimate in the light of the Clausius-Duhem inequality, a statementof the 2nd law of thermodynamics (e.g., Sec. 3.2 in Lubliner, 2008). Rates andamounts of strain softening in rocks are still to be better constrained but reducingthe uncertainty associated with them is an independent issue.In this paper, we first derive a set of governing equations for thermo-mechanicallycoupled tectonic systems. We start from the generic thermodynamic principlesclosely following the procedure of (Wright, 2002) to derive the governing equationthat incorporate the three elements discussed above. We focus on the elasto-visco-plastic material that has compressible elasticity, thermal stress and strain wakeningplasticity. Integration of such kinematic and constitutive models into the generalenergy balance will be crucial for realistically modeling geological systems. Wethen implement the derived governing equations in DES3D, an open source fi-nite element code for geodynamic modeling, and verify the implementation semi-analytically. Finally, we explore the e ff ects of thermo-mechanical coupling on thelong-term evolution of large-o ff set normal faults with emphasis on the role of vol-umetric inelastic strain. Since extensibly studied and well understood, the normalfault systems would allow us to isolate the new e ff ects introduced by the coupled4hysics.
2. Derivation of governing equations
Several theoretical works have derived a set of governing equations for a thermo-mechanical system from the general form of the thermodynamic principles. Thecommon procedure is to relate the rate of change of the internal energy appearingin the statement of energy balance to that of thermodynamic potentials such as theHelmholtz free energy and the Gibbs free energy. Thermodynamic potentials in-volve the capacity to do mechanical work in addition to heat content. For instance,the Helmholtz free energy is “the portion of the internal energy available for do-ing work at constant temperature” (p.263 in Malvern, 1969). The main di ff erencebetween the the Helmholtz and the Gibbs free energy is whether strain is an inde-pendent variable as in the former or stress is as in the latter. Since the definitionsof these thermodynamic potentials involve the product of temperature and entropy,the energy balance principle takes an intermediate form involving the time deriva-tive of entropy. The last step in deriving the temperature evolution equation is toexpress the time derivative of entropy in terms of that of temperature and othervariables. ? start from the energy conservation principle stated in terms of theHelmholtz free energy to derive the partial di ff erential equation for temperatureevolution as well as other equations that are coupled with it (e.g., mass conserva-tion equation and constitutive relations) for shear zone-developing systems in theEarth and other planets. Their final system of equations can consistently describethe feedback processes among energy, rheology and other variables such as grainsize and water content in shear zone formation. Similarly, Lyakhovsky et al. (1997)show how an evolution equation for elastic damage can be derived from the energy5onservation principle. Also starting from the energy balance equation in termsof the Helmholtz free energy, they include a non-dimensional variable quantifyingthe amount of damage along with temperature and infinitesimal elastic strain asindependent variables of the free energy. By treating damage process as a sourceof entropy in the intermediate equation for entropy evolution, they derive a dam-age evolution equation that is proportional to the rate of free energy change withdamage.For completeness, we derive a temperature evolution equation from the genericenergy balance principle involving the Gibbs free energy ( g ), following Wright(2002). This form of the energy balance principle is di ff erent from those of re-lated studies in geodynamics (Regenauer-Lieb and Yuen, 2003; Lyakhovsky et al.,1997) in that those earlier studies used another thermodynamic potential, the Helmholtzfree energy. In the context of continuum thermodynamics, the Gibbs free energyis a function of Cauchy stress ( σ ), temperature ( T ) and a set of internal variables( q j , j = , , ... n ) while the Helmholtz energy has elastic strain in place of Cauchystress. The work by Wright (2002) inspired our choice of thermodynamic potentialbut one di ff erence is that we assume infinitesimal strain while Wright (2002) usedthe finite strain kinematics.Additive strain rate decomposition is allowed under our assumption of in-finitesimal strain. Furthermore, the class of material we are interested in is elasto-visco-plastic. Under these assumptions, total strain rate is decomposed as follows:˙ ǫ = ˙ ǫ e + ˙ ǫ p + ˙ ǫ v , (1)where ˙ ǫ e , ˙ ǫ p , and ˙ ǫ v are elastic, plastic and viscous strain rate, respectively.Following Wright (2002), we define the Gibbs free energy per unit mass as g ( σ , T , q j ) = − e + T η + ρ σ : ǫ e , (2)6here e is the internal energy per unit mass, η is entropy per unit mass, ρ is densityand ǫ e is elastic strain tensor. Only the elastic strain appears in the definition be-cause it is the only component of strain that can contribute to stored energy. Takingthe totral di ff erential of g and rearranging terms, we get the di ff erential of internalenergy: de = − dg + η dT + T d η − ρ σ : ǫ e d ρ + ρ ǫ e : d σ + ρ σ : d ǫ e . (3)Since g = g ( σ , T , q j ), the total di ff erential of g is also dg = ∂ g ∂ σ : d σ + ∂ g ∂ T dT + n X j = ∂ g ∂ q j dq j , (4)where q = ρ according to the conventions of Wright (2002). Elastic strain andentropy (c.f., Wright, 2002) are defined as ǫ e ≡ ρ ∂ g ∂ σ and η ≡ ∂ g ∂ T . (5)From the equations (3), (4) and (5) we get de = − n X j = ∂ g ∂ q j dq j + T d η − ρ σ : ǫ e d ρ + ρ σ : d ǫ e (6)Since q = ρ , the terms containing dq and d ρ can be grouped together. With thegrouping, equation (6) becomes de = ρ σ : d ǫ e + T d η − ∂ g ∂ρ + ρ σ : ǫ e ! d ρ − n X j = ∂ g ∂ q j dq j . With two new notations Q ≡ − ρ ∂ g ∂ρ + ρ σ : ǫ e ! (7)and Q j ≡ − ρ ∂ g ∂ q j (j =
1, 2, 3, . . . ) , (8)7quation (6) is simplified to de = ρ σ : d ǫ e + T d η + ρ n X j = Q j dq j . (9)Di ff erentiating the equation (9) with respect to time ( t ) we have dedt = T d η dt + ρ σ : ˙ ǫ e + ρ n X j = Q j dq j dt . (10)Multiplying the equation (10) by ρ , we get the following equation for the time rateof change of internal energy per volume: ρ dedt = ρ T d η dt + σ : ˙ ǫ e + n X j = Q j dq j dt . (11)To relate the material time derivative of the internal energy given in (11) tothe energy balance principle, we recall the general form of the energy balanceequation (e.g., Kennett and Bunge, 2008; Malvern, 1969; Wright, 2002): ρ dedt = σ : ∇ v − ∇ · q + ρ s , (12)where s is a heat energy source or sink per mass, q is heat flux and v is velocity.According to the additive decomposition of strain rate in (1), σ : ˙ ǫ e = σ : ˙ ǫ − σ : (˙ ǫ p + ˙ ǫ v ) . (13)Since the double dot product of the symmetric σ and the anti-symmetric part of ∇ v is zero, σ : ∇ v = σ : ˙ ǫ . By eliminating the time derivative of internal energy fromequations (11) and (12) and then using (13), we get ρ T d η dt + σ : ˙ ǫ − σ : (˙ ǫ p + ˙ ǫ v ) + n X j = Q j dq j dt − σ : ˙ ǫ − ρ s + ∇ · q = , which is simplified to ρ T d η dt = σ : (˙ ǫ p + ˙ ǫ v ) − ∇ · q + ρ s − n X j = Q j dq j dt . (14)8he final step to derive a partial di ff erential equation for temperature from (14)is to relate the time derivative of entropy per mass to temperature. The equipres-ence principle (Malvern, 1969) requires the entropy η to have the same set ofindependent variables as the Gibbs free energy. In other words, the entropy perunit mass is also a function of the Cauchy stress ( σ ), the temperature ( T ) and a setof internal variables ( q j , j = , , ... n ). The total di ff erential of η = η ( σ , T , q j ) is d η = ∂η∂ σ : d σ + ∂η∂ T dT + n X j = ∂η∂ q j dq j . Identifying the first internal variable ( q ) with density ( ρ ) again, we get d η = ∂η∂ σ : d σ + ∂η∂ T dT + ∂η∂ρ d ρ + n X j = ∂η∂ q j dq j . (15)When di ff erentiated with respect to time, the equation (15) becomes d η dt = ∂η∂ σ : d σ dt + ∂η∂ T dTdt + ∂η∂ρ d ρ dt + n X j = ∂η∂ q j dq j dt . (16)According to (5), ∂η∂ σ = ∂∂ σ ∂ g ∂ T = ∂∂ T ∂ g ∂ σ = ρ ∂ ǫ e ∂ T . (17)From the definition of the specific heat at constant stress, c σ , we get the followingidentity: c σ ≡ ∂ e ∂ T = ∂∂ T − g + T η + ρ σ : ǫ e ! = T ∂η∂ T + ρ σ : ∂ ǫ e ∂ T . (18)Using (18), we get ∂η∂ T = T c σ − ρ σ : ∂ ǫ e ∂ T ! . (19)For convenience, we express partial derivatives of entropy per mass with respect todensity and other internal variables in terms of Q and Q j ( j = , , , . . . ), whichare defined in (7) and (8). ∂η∂ρ = ∂∂ρ ∂ g ∂ T = ∂∂ T ∂ g ∂ρ = − ρ ∂ Q ∂ T − ρ σ : ∂ ǫ e ∂ T . (20)9nd ∂η∂ q j = ∂∂ q j ∂ g ∂ T = ∂∂ T ∂ g ∂ q j = − ρ ∂ Q j ∂ T , (21)for j = , , , . . . . Plugging (17), (19), (20) and (21) and into the equation (16), weget d η dt = ρ ∂ ǫ e ∂ T : d σ dt + T c σ − ρ σ : ∂ ǫ e ∂ T ! dTdt + − ρ ∂ Q ∂ T − ρ σ : ∂ ǫ e ∂ T ! d ρ dt − ρ n X j = ∂ Q j ∂ T dq j dt . This equation can be further simplified to d η dt = T c σ − ρ σ : ∂ ǫ e ∂ T ! dTdt + ρ ∂ ǫ e ∂ T : d σ dt − ρ σ : ∂ ǫ e ∂ T ! d ρ dt − ρ n X j = ∂ Q j ∂ T dq j dt . (22)Substituting the equation (22) into the equation (14) we get: ρ T T c σ − ρ σ : ∂ ǫ e ∂ T ! dTdt + ρ ∂ ǫ e ∂ T : d σ dt − ρ σ d ρ dt ! − ρ n X j = ∂ Q j ∂ T dq j dt − σ : (˙ ǫ p + ˙ ǫ v ) + n X j = Q j dq j dt − ρ s + ∇ · q = , which is simplified to ρ c σ − ρ σ : ∂ ǫ e ∂ T ! dTdt + T ∂ ǫ e ∂ T : d σ dt − ρ σ d ρ dt ! + n X j = Q j − T ∂ Q j ∂ T ! dq j dt − σ : (˙ ǫ p + ˙ ǫ v ) − ρ s + ∇ · q = . (23)Strain due to temperature change is assumed to be isotropic: i.e., ∂ ǫ e /∂ T = α l I ,where I is the identity tensor and α l is linear thermal expansion coe ffi cient and 1 / ffi cient, α v . Under this assumption, the10quation (23) is rearranged as follows: ρ c σ dTdt − σ : α l I dTdt = − ∇ · q + ρ s + σ : (˙ ǫ p + ˙ ǫ v ) − T α l I : d σ dt − ρ d ρ dt σ ! − n X j = Q j − T ∂ Q j ∂ T ! dq j dt . (24)Using the identity σ : α l I = − α l p = − α v p where p = − σ ii , we have thefollowing final form of the energy balance equation:( ρ c p + p α v ) dTdt = − ∇ · q + ρ s + σ : (˙ ǫ p + ˙ ǫ v ) + T α v d pdt − p T α v ρ d ρ dt − n X j = Q j − T ∂ Q j ∂ T ! dq j dt . (25)The last term of (25) corresponds to the changes in energy due to internal vari-ables only. These terms are often parametrized into a coe ffi cient for the inelasticpower term as in χ σ : (˙ ǫ p + ˙ ǫ v ) (Regenauer-Lieb and Yuen, 2003; Wright, 2002).With this parametrization, the equation (25) is simplied to( ρ c p + p α v ) dTdt = −∇ · q + ρ s + χ σ : (˙ ǫ p + ˙ ǫ v ) + T α v d pdt − p T α v ρ d ρ dt . (26) χ close to 1 means that the energy change due to changes in internal variables isnegligible relative to inelastic power and most of the inelastic power is convertedto heat production. To our knowledge, the value of χ is not very well constrainedfor rocks but at least for metals, it is either close to 1 or saturates towards 1 withincreasing plastic strain (Wright, 2002). Regenauer-Lieb and Yuen (2003) used0.85 as the value of χ . χ is assumed to be 1 in this study. The form of the mass conservation equation involving the material derivativeis d ρ dt = − ρ ∇ · v . (27)11n the Lagrangian description of motion that we are going to adopt, the above timederivative is understood as the partial derivative with respect to time, not as thematerial time derivative. We substitute (27) into the last term of (26) to get( ρ c p + p α v ) dTdt = −∇ · q + ρ s + σ : (˙ ǫ p + ˙ ǫ v ) + T α v d pdt + p T α v ∇ · v . (28) Assuming the linear isotropic elasticity and temperature change δ T , we use thefollowing thermoelastic constitutive equations (e.g., Boley and Weiner, 1997): σ i j = λε kk δ i j + G ε i j − K α v δ T , (29)where λ and G are the Lam´e’s constant and K is the bulk modulus defined as λ + G /
3. These thermal elastic constitutive equations become the basis for vis-coelastic or elastoplastic constitutive models as described in (Choi et al., 2013b).
Using a non-associated, strain-weakening plasticity for inducing shear bands asin (Choi et al., 2013b) has been described as if it violated the 2nd law of thermody-namics (e.g., Regenauer-Lieb and Yuen, 2003; Regenauer-Lieb et al., 2006, 2008).An accompanying criticism on the strain-weakening plasticity is that the softeningrules often used in tectonic modeling lack reliable constraints. While agreeing thatthe lack of experimental and observational constraints on the strain-softening rulesis problematic, we would like to point out that this problem does not invalidatethe approach itself. In fact, shear band formation induced by strain softening isperfectly legitimate in the light of the 2nd law of thermodynamics. The source ofprevious confusion must be in the erroneous notion that only a strain hardeningplasticity is “stable” in the sense that the plastic dissipation is positive therefore the12nd law of thermodyanmics is satisfied. However, neither strain softening nor non-associated flow rule necessarily violates the maximum plastic dissipation postulateas shown by Lubliner (2008, see Sec. 3.2.)The maximum plastic dissipation postulate requires that the yield surface beconvex in the principal stress space and the isotropic strain-softening plasticitymodel often used in tectonic modeling! satisfy this requirement by maintainingthe convexity. For instance, when a strain-softening in the Mohr-Coulomb plasticmodel is realized through reduction in cohesion and friction angle, the hexagonalcone-shaped, therefore convex, yield surface in the principal stress space alwaysremain convex except in the degenerate case where both friction angle and cohesionare zero.Under a typical setting for numerical tectonic modeling in which plastic flowis incompressible, i.e., dilation angle is 0 ◦ , and the friction angle is about 30 ◦ , theprincipal stresses and plastic strain rates are indeed “non-coaxial” (Regenauer-Lieb and Yuen,2003; Regenauer-Lieb et al., 2006). However, this fact does not necessarily meanthat such a non-associated flow rule violates the maximum plastic dissipation. Fric-tion angles usually assumed in tectonic models are never much greater than 30 ◦ anddilation angles are less than that. Since principal stress and plastic strain orienta-tions are normal to a yield surface and a flow potential, respectively, the di ff erencebetween their orientations is equal to that of the friction and the dilation angle.Because the double contraction of stress and plastic strain rate in the definitionof plastic dissipation rate is equivalent to the projection operation, an acute anglebetween the principal stresses and plastic strain rates guarantees a positive plas-tic dissipation that indicates the satisfaction of the maximum plastic dissipationpostulate.These considerations validate our adoption of the strain-weakening, non-associated,Mohr-Coulomb plastic model. 13 . Benchmarks We verify the governing equations implemented into
DES3D , an unstructuredfinite element code for long-term tectonic modeling (Choi et al., 2013b), usingthree benchmark problems. The benchmarks are derived from the standard oe-dometer test (e.g., Davis and Selvadurai, 2002; Choi et al., 2013b). A cubic blockof the Mohr-Coulomb plastic material is compressed in one direction while mo-tions in the other directions are restricted (Fig. 1). The symmetry and boundaryconditions of the problem make it su ffi cient to discretize a square domain withtwo triangular elements. The simplicity of the problem allows for at least semi-analytic solutions. In benchmark-1, we include only the plastic power ( σ : ˙ ǫ p ) asa heat source and ignore the di ff usion term. The density is assumed to be constant.Benchmark-2 solves the following equation( ρ c p + p α v ) dTdt = σ : ˙ ǫ p + T α v d pdt + p T α v ∇ · v , (30)and the density is updated according to eq. (27). Benchmark-3 verifies the thermalstress calculations in DES3D for a uniform temperature that increases linearly intime. We intentionally use a low density (1 kg / m ) to make the contribution fromplastic power non-negligible. Parameters used in the benchmarks are listed in Table1. 14 = 1 m L = m v x = - m / s Free slip
Figure 1: Model geometry and boundary conditions for the plastic oedometer test. This figure andassociated running / plotting scripts available under Ahamed and Choi (2016) The model geometry and boundary conditions give the following componentsof total strain rate (˙ ε ): ˙ ε xx = v x L + v x t , (31)˙ ε yy = , (32)˙ ε zz = , (33)where t is time, v x is the boundary velocity set to be − − m / s and L is the edgelength of the cube equal to 1 m. 15 able 1: Model parameters for Benchmarks Paramter Symbol ValueBulk Modulus K
200 MPaShear Modulus G
200 MPaCohesion C φ ◦ Dilationa Angle Ψ ◦ Initail Temperature T
273 KReference density ρ / m Volumetric expansion coe ffi cient α − The components of total strain ( ε ) are given as ε xx = ln (cid:18) L + v x tL (cid:19) , (34) ε yy = , (35) ε zz = . (36) Before yielding, stresses are updated according to the linear isotropic elasticity.In terms of the Lam´e’s constants ( λ , G ), stress components are given as: σ xx = ( λ + G ) ln (cid:18) L + v x tL (cid:19) , (37) σ yy = σ zz = λ ln (cid:18) L + v x tL (cid:19) . (38)16 .1.3. Mohr-Coulomb yield function and flow potential We use the following forms of the Mohr-Coulomb yield function ( f ) and flowpotential ( g ): f ( σ xx , σ yy ) = σ xx − N φ σ yy + C p N φ , (39) g ( σ xx , σ yy ) = σ xx − + sin ψ − sin ψ σ yy . (40)where N φ is defined as (1 + sin φ ) / (1 − sin φ ), φ is the firction angle, ψ is the dilationangle and C is the cohesion. cr ) We denote the time when yielding occurs for the first time as t cr . Since f ( σ xx , σ yy ) = t = t cr , ( λ + G ) ln (cid:18) L + v x t cr L (cid:19) − N φ λ ln (cid:18) L + v x t cr L (cid:19) + C p N φ = . (41)Solving the above equation for t cr , we get t cr = Lv x exp − C p N φ ( λ + G ) − N φ λ − . (42) We get the following expressions for plastic strain for t ≥ t cr : ε pxx = β ( t ) , (43) ε pyy = − N ψ β ( t ) , (44) ε pzz = − N ψ β ( t ) , (45)where β ( t ) is the plastic consistency paramter to be determined. As in the classicalplasticity theory (e.g., Lubliner, 2008), the consistency parameter is zero beforeyielding ( t ≤ t cr ). 17he plastic strain rates are similarly defined in terms of the rate of consistencyparameter ( ˙ β ( t )): ˙ ε pxx = β ( t ) , (46)˙ ε pyy = − N ψ ˙ β ( t ) , (47)˙ ε pzz = − N ψ ˙ β ( t ) . (48) The above definitions of total and plastic strain lead to the following expres-sions for elastic strain as a function of time: ε exx = ε xx − ε pxx = ln (cid:18) L + v x tL (cid:19) − β ( t ) , (49) ε eyy = ε yy − ε pyy = N ψ β ( t ) , (50) ε ezz = N ψ β ( t ) . (51) β ( t )After yielding occurs, stresses are updated as follows: σ xx = ( λ + G ) ε exx + λ ( ε eyy + ε ezz ) , (52) σ yy = ( λ + G ) ε eyy + λ ( ε ezz + ε exx ) , (53) σ zz = ( λ + G ) ε ezz + λ ( ε exx + ε eyy ) . (54)(55)18lugging (49) to (51) into the above equations, we get σ xx ( t ) = ( λ + G ) (cid:18) ln (cid:18) L + v x tL (cid:19) − β ( t ) (cid:19) + λ N ψ β ( t ) = ( λ + G ) ln (cid:18) L + v x tL (cid:19) − h λ + G ) − λ N ψ i β ( t ) , (56) σ yy ( t ) = ( λ + G ) N ψ β ( t ) + λ (cid:18) N ψ β ( t ) + ln (cid:18) L + v x tL (cid:19) − β ( t ) (cid:19) = λ ln (cid:18) L + v x tL (cid:19) + h λ + G ) N ψ − λ i β ( t ) , (57) σ zz ( t ) = σ yy ( t ) . (58)For t ≥ t cr , these stresses should always satisfy the yield condition: σ xx − N φ σ yy + C p N φ = . (59)Substituting (56) and (57) for the stress components in the yield condition, we get( λ + G ) ln (cid:18) L + v x tL (cid:19) − h λ + G ) − λ N ψ i β ( t ) − N φ λ ln (cid:18) L + v x tL (cid:19) − h λ + G ) N φ N ψ − λ N φ i β ( t ) + C p N φ = . (60)Solving the above equation for β ( t ), we get β ( t ) = h ( λ + G ) − N φ λ i ln (cid:16) L + v x tL (cid:17) + C p N φ λ + G ) N φ N ψ + λ + G ) − λ ( N φ + N ψ ) . (61)The time derivative of β ( t ) is given as˙ β ( t ) = ( λ + G ) − N φ λ λ + G ) N φ N ψ + λ + G ) − λ ( N φ + N ψ ) v x L + v x t . (62)19 .1.8. Pressure and pressure rate Since p = − ( σ xx + σ yy + σ zz ) /
3, we get the following from (56), (57) and (58): p ( t ) = − K ln (cid:18) L + v x tL (cid:19) + λ + G ! (1 − N ψ ) β ( t ) = − K ln (cid:18) L + v x tL (cid:19) + K (1 − N ψ ) β ( t ) , ˙ p ( t ) = − K v x L + v x t + λ + G ! (1 − N ψ ) ˙ β ( t ) = − K v x L + v x t + K (1 − N ψ ) ˙ β ( t ) , (63)where K is the bulk modulus defined as λ + / G . ∇ · v = ˙ ε xx + ˙ ε yy + ˙ ε zz = v x L + v x t . (64) Using the expression for ∇ · v , the mass balance equation becomes d ρ dt = − ρ ∇ · v = − v x L + v x t ! ρ. (65)Solving for ρ ( t ), we get ρ ( t ) = ρ LL + v x t , (66)where ρ is the reference density at t = The equation we are going to solve is ρ c p dTdt = σ : ˙ ε p , (67)20 .0 0.5 1.0 1.5 2.0 Displacement (x10 −2 m) T e m p e r a t u r e ( K ) AnalyticDES3D
Displacement (x10 −2 m) S t r e ss ( M P a ) |σ xx | Analytic DES3D|σ yy | Analytic DES3D Figure 2: Temperature ( T ), Stress xx ( | σ xx | ), and Stress yy ( | σ yy | ) are plotted against x-displacementfrom analytic and numerical solutions by DES3D for benchmark-1. This figure and associated run-ning / plotting scripts available under Ahamed and Choi (2016) where σ : ˙ ε p = σ xx ˙ ε pxx + σ yy ˙ ε pyy + σ zz ˙ ε pzz = β ( t ) σ xx ( t ) − N ψ ˙ β ( t ) σ yy ( t ) = β ( t ) h σ xx ( t ) − N ψ σ yy ( t ) i . (68)Since ρ and c p are assumed to be constant, dTdt = ρ c p σ : ˙ ε p = β ( t ) ρ c p h σ xx ( t ) − N ψ σ yy ( t ) i . (69)We use the forward Euler scheme to integrate the above time derivative of T withthe initial condition T (0) =
273 K.Results of the test are shown in Fig. 2. The relative error of the DES3D solutionrelative to the semi-analytic one is 0.004 %.21 .0 0.5 1.0 1.5 2.0
Displacement (x10 −2 m) T e m p e r a t u r e ( K ) AnalyticDES3D
Displacement (x10 −2 m) D e n s i t y ( k g / m ) Figure 3: Temperature ( T ) and Density ( ρ ) are plotted against x-displacement from analytic andnumerical solutions by DES3D for benchmark-2. This figure and associated running / plotting scriptsavailable under Ahamed and Choi (2016) The equation to solve is( ρ c p + p α v ) dTdt = σ : ˙ ε p + T α v d pdt + p T α v ∇ · v . (70)Using the derived expressions for the involved quantites and the forward Eulerscheme, we can integrate the above equation for the initial condition, T (0) = . We verify the implementation of the thermoelastic constitutive equations inDES3D. For simplicity, we prescribe a uniform temperature field, which is initially2273 K and increases linearly in time as T ( t ) = b t , (71)where b is a constant. We set b to be 0.4 K / s.Under this setting, we can analytically derive the elastic and plastic stress so-lutions for the oedometer test involving thermal stress. In the elastic regime, thestresses are given as σ xx = ( λ + G ) ε exx + λ ( ε eyy + ε ezz ) − K α v b t , (72) σ yy = ( λ + G ) ε eyy + λ ( ε ezz + ε exx ) − K α v b t , (73) σ zz = ( λ + G ) ε ezz + λ ( ε exx + ε eyy ) − K α v b t . (74)Using the expressions, (56), (57) and (58), we get σ xx ( t ) = ( λ + G ) (cid:18) ln (cid:18) L + v x tL (cid:19) − β ( t ) (cid:19) + λ N ψ β ( t ) − K α v b t (75) = ( λ + G ) ln (cid:18) L + v x tL (cid:19) − h λ + G ) − λ N ψ i β ( t ) − K α v b t , (76) σ yy ( t ) = ( λ + G ) N ψ β ( t ) + λ (cid:18) N ψ β ( t ) + ln (cid:18) L + v x tL (cid:19) − β ( t ) (cid:19) − K α v b t (77) = λ ln (cid:18) L + v x tL (cid:19) + h λ + G ) N ψ − λ i β ( t ) − K α v b t , (78) σ zz ( t ) = σ yy ( t ) . (79)We follow the procedure in Sec. 3.1.4 to compute the timing of the first yielding, t cr ,but use the Newton method because the equation for t cr does not allow a solution interms of elementary functions. Following the procedure for determining β ( t ) and˙ β ( t ) described in Sec. 3.1.7, we get β ( t ) = h ( λ + G ) − N φ λ i ln (cid:16) L + v x tL (cid:17) + ( N φ − K α v b t + C p N φ λ + G ) N φ N ψ + λ + G ) − λ ( N φ + N ψ ) , (80)and ˙ β ( t ) = h ( λ + G ) − N φ λ i v x L + v x t + ( N φ − K α v b λ + G ) N φ N ψ + λ + G ) − λ ( N φ + N ψ ) . (81)23s before, β ( t ) = ˙ β ( t ) = t < t cr .Di ff erential stress arising in this benchmark is the same with those of theisothermal case but pressure in this test is greater due to the compressional ther-mal stress caused by the prescribed temperature increase. As a result, the firstyielding in benchmark-3 should occur at a greater value of di ff erential stress, orequivalently, displacement than in the isothermal case. In other words, than in theisothermal case. The results of benchmark-3 shows in Fig. 4 are consistent withthis expectation. The relative errors of σ xx and σ yy computed with DES3D are0.0028 % and 0.01 %.
4. Discussion
Using the verified implementation of the “full” energy balance equation, (28),as well as the thermo-elasticity, we now systematically assess the impact of thethermo-mechanically coupled governing equations on a more geologically rele-vant problem. In this study, we choose the large-o ff set normal fault evolution (e.g.Lavier et al., 2000) as an example. Studied extensively and understood well, thissystem will facilitate identification and attribution of di ff erences between the mod-els of the newly-introduced physics and the uncoupled, isothermal ones.We create a normal fault in an extending Mohr-Coulomb elasto-plastic layerwhich is initially 100 km long and 10 km thick (Fig. 5). Both sides of the layerare pulled at a constant velocity of 0.5 cm / yr. The bottom boundary is supportedby the Winkler foundation (Watts, 2001, pp.95). To induce strain localization, wedecrease cohesion from 40 MPa to 8 MPa linearly as plastic strain increases to 1.We add a weak zone at the bottom center of the domain to initiate a fault fromthere. We impose topographic smoothing of the di ff usion type with a transportcoe ffi cient of 10 − m / s (Turcotte and Schubert, 2014, pp. 225).24 .0 0.5 1.0 1.5 2.0 Displacement (x10 −2 m) S t r e ss ( M P a ) |σ xx | Analytic DES3D Isothermal|σ yy | Analytic DES3D Isothermal Figure 4: | σ xx | and | σ yy | from the analytic (circles) and the DES3D numerical solution (solid lines) forbenchmark-3, plotted against displacement. The corresponding analytic solutions for the isothermalcase (dashed lines) are shown for comparison. This figure and associated running / plotting scriptsavailable under Ahamed and Choi (2016) inkler Foundation k m
100 km . c m / y ea r . c m / y ea r Density = 2700 Kg/m Weak inhomogeneity
Figure 5: Setup of the large o ff set normal fault models. This figure and associated running / plottingscripts available under Ahamed and Choi (2016) We set up five models that di ff er in the form of the energy balance equationand in the presence of volumetric plastic strain. In model-1, we do not considerany heat-generating mechanism. As a result, the initial temperature, which is uni-formly 0 ◦ C, does not change in time. The behavior of model-1 should be similar tothat of the “unlimited, footwall-snapping” mode of (Lavier et al., 2000). In the restof models (model-2 to 5) we solve the full energy balance equation (28) with theheat source / sink term, s , equal to zero and with the temperature feedback to rheol-ogy through thermo-elasticity. To explore the e ff ects of volumetric plastic strain,we use a non-zero dilation angle in model-3 to 5. The dilation angle is fixed inmodel-3 while we reduce it to zero as the accumulated plastic strain increases to aprescribed value of 1.0 in model-4 and 5. The model-5 is the same as model-4 butwe intentionally include only the deviatoric part of the plastic power ( σ : ˙ ε p ) inequation (28). By comparing model-4 and 5, we can decide whether the volumet-ric plastic power can be ignored. Table 2 summarizes the di ff erences among themodels. Table 3 shows the list of the parameters used in the models.We see noticeable di ff erences in fault geometry and shape of the core complexbetween model-1 and model-2 that solves the full energy balance equation andincludes thermal stresses. The overall behavior of the faults in model-1 is similarto the unlimited, footwall snapping mode of Lavier et al. (2000). The geometry26 able 2: Description of the normal fault evolution models Models Energy balance Dilation angle Plastic powerModel-1 Heat di ff usion only 0 ◦ N / AModel-2 Full 0 ◦ TotalModel-3 Full 10 ◦ TotalModel-4 Full 10 ◦ to 0 ◦ TotalModel-5 Full 10 ◦ to 0 ◦ Deviatoric
Table 3: Parameters for the normal fault evolution models
Parameter Symbol ValueBulk Modulus K
50 GPaLam´e’s constant λ
30 GPaShear Modulus G
30 GPaInitial Cohesion C
40 MPaFriction Angle φ ◦ Dilation Angle Ψ ◦ Density ρ Kg / m Volumetric expansion coe ffi cient α K − of the primary fault of model-2 is almost the same as that of the model-1 untilabout 17 km of extension. Around this time, however, model-2 forms a secondaryfault but model-1 does not yet. After 20 km of extension, more secondary faultsstart form in both models but their geometry and location are not identical. Forinstance, a fault block forms in model-2 after about 30 km of extension when anew high-angle fault forms next to the first secondary fault. Model-1 does notdevelop a corresponding structure including a fault block after the same amountof extension. Model-1 forms a single secondary fault initially, which eventually27onnects with the second secondary fault after about 40 km of extension. At 30 kmof extension model-2 has three well-developed secondary faults while model-1 hastwo. Initially, the shape of core complex remains same until 10 km of extension.Because of di ff erent fault geometry the shape of the core complex started to di ff erfrom each other from 20 km of extension. At 40 km of extension the core-complexof the model-2 becomes elongated with the almost flat surface while the model-1has more rounded one.Figure 7 shows the distribution of temperature in model-2 and the correspond-ing thermal stresses after 10 and 30 km of extension. The temperature increasenear the active faults are greater than 15 ◦ C and the magnitude of thermal stressis as large as 40 MPa. The thermal stress magnitude is a non-negligible fraction,about 10 %, of lithostatic stress near the bottom. Thus, the di ff erent behaviors ofthe models can be attributed to the presence of thermal stresses only in model-2. a Model-1 b Model-2 Plastic strain0 1010 Km20 Km30 Km40 Km Primary fault Primary faultSecondary faultSecondary fault Fault blockFault block
Figure 6: Plastic strain distribution of (a) model-1 and (b) model-2 after 10, 20, 30 and 40 km ofextensions. This figure and associated running / plotting scripts available under Ahamed and Choi(2016) Identical with model-2 except a non-zero dilation angle of 10 ◦ , model-3 devel-ops an unrealistically elevated core-complex with a relief greater than 10 km andshear bands that are much wider than those in model-1 and 2. Expansion of shear28 o C) 25030 km of extension ba
30 km of extension
Figure 7: (a) Thermal stress and (b) temperature distribution of model-2 after 10 and 30 km ofextension. This figure and associated running / plotting scripts available under Ahamed and Choi(2016) bands when dilation angle is non-zero is kinematically expected. The dilation an-gle fixed at 10 ◦ sustains a higher level of plastic power than the non-dilational caseszero dilation angle. The greater plastic power leads to a greater amount of temper-ature change and thermal stresses. The thermal stresses pushes up the free-tractiontop surface, creating the highly elevated core-complex. Model-310 Km20 Km30 Km40 Km
Plastic strain0 10
Figure 8: Same with Fig. 6 but for model-3. This figure and associated running / plotting scriptsavailable under Ahamed and Choi (2016) Model-4 b Model-5
Plastic strain0 10
10 Km20 Km30 Km40 Km
Figure 9: Same with Fig. 6 but for (a) model-4 and (b) model-5. This figure and associated run-ning / plotting scripts available under Ahamed and Choi (2016) To isolate the e ff ect of volumetric plastic power, we construct the model-4 and5 that are identical except that model-4 considers the total plastic power whilemodel-5 includes only the deviatoric part of the plastic power. In order to preventthe unrealistic expansion of shear bands as in model-3, the dilation angle is reducedwith accumulated plastic strain in both model-4 and 5 (Choi and Petersen, 2015).Models-4 and 5 do not show notable di ff erences until 10 km extension (Figure9). However, after about 20 km of extension, they start to exhibit di ff erences in thefault geometry as well as in the timing of rider block formation. A rider block startsto form at 17 km of extension in model-4 but model-5 forms one only after morethan 37 km of extension. After 40 km of extension, model-5 has three secondaryfaults produced by the footwall snapping but model-4 still has only two secondaryfault. The topographic relief in model-4 and 5 is about 5 km, which is much greaterthan the relief of the previous models for the large o ff set normal fault (Lavier et al.,2000; Choi et al., 2013a) but about half of that of model-3 with a fixed dilationangle of 10 ◦ . The excessive reliefs in the models with non-zero dilation angle,model-3 to 5, suggest that the rate of dilation angle reduction should be greaterthan that of this study to have 1-2 km topographic relief as in the previous studies30n large-o ff set normal faults. The noticeable di ff erences between model-4 and 5suggest to include the total plastic power in the energy balance rather than onlythe deviatoric component. Ignoring the volumetric plastic power as in model-5 isinconsistent with the kinematics in the first place.
5. Conclusions
We derive a temperature evolution equation based on the energy balance prin-ciple that accounts for deformation-related energy changes as well as the conven-tional heat energy di ff usion and advection. An explicit-time finite element solu-tion procedure for the derived equation is implemented in DES3D, an unstructuredmesh finite element solver for geodynamics. We verify the implementation usingthree benchmark tests that have semi-analytic solutions. Benchmark-1 includesonly the total plastic power term in the energy balance equation. Benchmark-2 in-cludes the full energy equation as well as the mass balance equation. Benchmark-3prescribes a temperature change linear in time and computes thermal stresses. Thenumerical solutions from DES3D for all the benchmarks show an excellent matchwith the corresponding semi-analytic solutions. Coupling the verified implemen-tation of the full temperature evolution with the strain weakening Mohr-Coulombplastic rheology through thermal stresses, we also explore the long-term behaviorsof a large-o ff set normal fault. We find that the temperature arising mostly fromthe inelastic power term in the full energy balance has non-negligible e ff ects onthe long-term evolution of normal fault. For instance, when the plastic power isconsidered, temperature increases by more than 15 ◦ C and the magnitude of theassociated thermal stress is as large as 40 MPa. These extra stresses appear topromote the formation of secondary faults and an elongated core complex. Whenthe dilational plastic deformation is enabled, our models develop even more dif-31erences in faulting behaviors such as the formation of a rider block and the greattopographic relief, compared to the previous models for the large-o ff set normalfault with the same parameters and geometry. We conclude that the new coupledgoverning equations presented in this study has significant impact on long-termtectonics. Although not so strong as to cause a fundamental shift in faulting modesunder tested conditions, the e ff ects of deformation energetics might be criticallyimportant in other geologically relevant systems and conditions. ReferencesReferences
Ahamed, S., Choi, E., 2016. The e ff ects of deforma-tion energetics in Long-term tectonics modeling URL: https://figshare.com/articles/Figures/4213599 ,doi: .Biot, M.A., 1961. Theory of folding of stratified viscoelasticmedia and its implications in tectonics and orogenesis. Bul-letin of the Geological Society of America 72, 1595–1620.doi: .Boley, B.A., Weiner, J.H., 1997. Theory of thermal stresses. Dover Publications,Inc., Mineola, New York.Borja, B.R.I., Regueiro, R.A., Lai, T.Y., 2000. Fe Modeling of Strain Localizationin Soft Rock. Journal of Geotechnical and Geoenvironmental Engineering 126,335–343. doi: .Brace, W., Kohlstedt, D., 1980. Limits on lithospheric stress imposed by laboratoryexperiments. Journal of Geophysical Research: Solid Earth 85, 6248–6252.32race, W., Paulding, B., Scholz, C., 1966. Dilatancy in the fracture of crystallinerocks. Journal of Geophysical Research 71, 3939–3953.Choi, E., Buck, W.R., Lavier, L.L., Petersen, K.D., 2013a. Using corecomplex geometry to constrain fault strength. Geophysical Research Let-ters 40, 3863–3867. URL: http://doi.wiley.com/10.1002/grl.50732 ,doi: .Choi, E., Gurnis, M., 2008. Thermally induced brittle deformation in oceaniclithosphere and the spacing of fracture zones. Earth Planet. Sci. Lett. 269, 259–270. doi: .Choi, E., Petersen, K.D., 2015. Making coulomb angle-oriented shear bands innumerical tectonic models. Tectonophysics 657, 94–101.Choi, E., Tan, E., Lavier, L., Calo, V., 2013b. Dynearthsol2d: An e ffi cient unstruc-tured finite element method to study long-term tectonic deformation. Journal ofGeophysical Research: Solid Earth 118, 2429–2444.Choi, E.s., Lavier, L., Gurnis, M., 2008. Thermomechanics of mid-ocean ridgesegmentation. Physics of the Earth and Planetary Interiors 171, 374–386.Connolly, J., 2009. The geodynamic equation of state: what and how. Geochem-istry, Geophysics, Geosystems 10.Davis, R.O., Selvadurai, A.P.S., 2002. Plasticity and Geomechanics. CambridgeUniversity Press, Cambridge, UK, Cambridge, UK.Gerya, T.V., Yuen, D.A., Maresch, W.V., 2004. Thermomechanical modelling ofslab detachment. Earth and Planetary Science Letters 226, 101–116.33oetze, C., Evans, B., 1979. Stress and temperature in the bending lithosphere asconstrained by experimental rock mechanics. Geophysical Journal International59, 463–478.Het´enyi, G., Godard, V., Cattin, R., Connolly, J.A., 2011. Incorporating meta-morphism in geodynamic models: the mass conservation problem. GeophysicalJournal International 186, 6–10.Hobbs, B.E., Ord, A., Regenauer-Lieb, K., 2011. The thermodynamics of de-formed metamorphic rocks: a review. Journal of Structural Geology 33, 758–818.Hunsche, U., 1991. Volume change and energy dissipation in rock salt duringtriaxial failure tests, in: Mechanics of Creep Brittle Materials 2. Springer, pp.172–182.Hunt, G.W., Wadee, M.K., 1991. Comparative Lagrangian Formulationsfor Localized Buckling. Proceedings of the Royal Society A: Math-ematical, Physical and Engineering Sciences 434, 485–502. URL: http://rspa.royalsocietypublishing.org/cgi/doi/10.1098/rspa.1991.0109 ,doi: .Hyndman, R.D., Peacock, S.M., 2003. Serpentinization of the forearc mantle.Earth and Planetary Science Letters 212, 417–432.Ita, J., King, S.D., 1994. Sensitivity of convection with an endothermic phasechange to the form of governing equations, initial conditions, boundary condi-tions, and equation of state. Journal of Geophysical Research: Solid Earth 99,15919–15938. 34ennett, B., Bunge, H., 2008. Geophysical continua. Cambridge University, Kap, 10–11.Korenaga, J., 2007. Thermal cracking and the deep hydration of oceanic litho-sphere: A key to the generation of plate tectonics? Journal of GeophysicalResearch: Solid Earth 112, 1–20. doi: .Lavier, L.L., Buck, W.R., Poliakov, A.N., 2000. Factors controlling normal faulto ff set in an ideal brittle layer. Journal of Geophysical Research: Solid Earth(1978–2012) 105, 23431–23442.Lubliner, J., 2008. Plasticity Theory. Dover Publications, Inc., Mineola, New York.Lyakhovsky, V., Ben-Zion, Y., Agnon, A., 1997. Distributed damage, faulting, andfriction. Journal of Geophysical Research: Solid Earth 102, 27635–27649.Malvern, L.E., 1969. Introduction to the Mechanics of a Continuous Medium.Prentice-Hall, Inc., Upper Saddle River, New Jersey.Ranalli, G., Murphy, D.C., 1987. Rheological stratification of the lithosphere.Tectonophysics 132, 281–295.Read, H., Hegemier, G., 1984. Strain softening of rock, soil and con-crete a review article. Mechanics of Materials 3, 271–294. URL: http://linkinghub.elsevier.com/retrieve/pii/0167663684900280 ,doi: .Regenauer-Lieb, K., Rosenbaum, G., Weinberg, R.F., 2008. Strain localisation andweakening of the lithosphere during extension. Tectonophysics 458, 96–104.Regenauer-Lieb, K., Weinberg, R.F., Rosenbaum, G., 2006. The e ff ect of energyfeedbacks on continental strength. Nature 442, 67–70.35egenauer-Lieb, K., Yuen, D., 2003. Modeling shear zones in geo-logical and planetary sciences: solid- and fluid-thermal–mechanicalapproaches. Earth-Science Reviews 63, 295–349. URL: http://dx.doi.org/10.1016/S0012-8252(03)00038-2 ,doi: .Turcotte, D.L., Schubert, G., 2014. Geodynamics. Cambridge University Press.Watts, A.B., 2001. Isostasy and Flexure of the Lithosphere. Cambridge UniversityPress.Wright, T.W., 2002. The physics and mathematics of adiabatic shear bands. Cam-bridge University Press.Yuen, D.A., Quareni, F., Hong, H.J., 1987. E ff ects from equation of state and rhe-ology in dissipative heating in compressible mantle convection. Nature 326,67–69. URL: ,doi:10.1038/326067a0