Incorporating Deformation Energetics in Long-Term Tectonic Modeling

TL;DR

This paper develops a comprehensive energy balance equation for long-term tectonic modeling, incorporating heat and mechanical work, and demonstrates its implementation and effects on fault evolution in numerical simulations.

Contribution

It introduces a full energy balance equation accounting for heat and mechanical work in tectonic models, implemented in DES3D, and explores its impact on fault system evolution.

Findings

Models with full energy balance produce more secondary faults.

Full energy considerations lead to elongated core complexes.

Persistent inelastic deformation significantly affects fault evolution.

Abstract

The deformation-related energy budget is usually considered in the simplest form or even completely omitted from the energy balance equation. We derive a full energy balance equation that accounts not only for heat energy but also for mechanical (elastic, plastic and viscous) work. The derived equation is implemented in DES3D, an unstructured finite element solver for long-term tectonic deformation. We verify the implementation by comparing numerical solutions to the corresponding semi-analytic solutions in three benchmarks extended from the classical oedometer test. Two of the benchmarks are designed to evaluate the temperature change in a Mohr-Coulomb elasto-plastic square governed by a simplified equation involving plastic power only and by the full temperature evolution equation, respectively. The third benchmark differs in that it computes thermal stresses associated with a prescribed uniform temperature increase. All the solutions from DES3D show relative error less than 0.1%. We also investigate the long-term effects of deformation energetics on the evolution of large offset normal faults. We find that the models considering the full energy balance equation tend to produce more secondary faults and an elongated core complex. Our results for the normal fault system confirm that persistent inelastic deformation has a significant impact on the long-term evolution of faults, motivating further exploration of the role of the full energy balance equation in other geodynamic systems.

Figures

Tables

Table 1. Table 1: Model parameters for Benchmarks

Paramter	Symbol	Value
Bulk Modulus	$K$	200 MPa
Shear Modulus	$G$	200 MPa
Cohesion	$C$	1 MPa
Friction Angle	$\varphi$	10∘
Dilationa Angle	$\Psi$	10∘
Initail Temperature	$T_0$	273 K
Reference density	${\rho }_0$	1 kg/m3
Volumetric expansion coefficient	$\alpha$	3.5 K-1

Table 2. Table 2: Description of the normal fault evolution models

Models	Energy balance	Dilation angle	Plastic power
Model-1	Heat diffusion only	0∘	N/A
Model-2	Full	0∘	Total
Model-3	Full	10∘	Total
Model-4	Full	10∘ to 0∘	Total
Model-5	Full	10∘ to 0∘	Deviatoric

Table 3. Table 3: Parameters for the normal fault evolution models

Parameter	Symbol	Value
Bulk Modulus	$K$	50 GPa
Lamé’s constant	$\lambda$	30 GPa
Shear Modulus	$G$	30 GPa
Initial Cohesion	$C$	40 MPa
Friction Angle	$\varphi$	${30}^{\circ }$
Dilation Angle	$\Psi$	${10}^{\circ }$
Density	$\rho$	2700 $Kg/m^3$
Volumetric expansion coefficient	$\alpha$	3.5 $K^{-1}$

Equations

$$\dot{\mathbf{ϵ}}={\dot{\mathbf{ϵ}}}_e+{\dot{\mathbf{ϵ}}}_p+{\dot{\mathbf{ϵ}}}_v,$$

$$g(\mathbf{\sigma },T,q_j)=-e+T\eta +\frac{1}{\rho }\mathbf{\sigma }:{\mathbf{ϵ}}_e,$$

$$de=-dg+\eta dT+Td\eta -\frac{1}{{\rho }^2}\mathbf{\sigma }:{\mathbf{ϵ}}_{e}d\rho +\frac{1}{\rho }{\mathbf{ϵ}}_e:d\mathbf{\sigma }+\frac{1}{\rho }\mathbf{\sigma }:d{\mathbf{ϵ}}_e.$$

$$dg=\frac{\partial g}{\partial \mathbf{\sigma }}:d\mathbf{\sigma }+\frac{\partial g}{\partial T}dT+\sum _{j=0}^n\frac{\partial g}{\partial q_j}d{q}_j,$$

$${\mathbf{ϵ}}_e\equiv \rho \frac{\partial g}{\partial \mathbf{\sigma }}\text{and}\eta \equiv \frac{\partial g}{\partial T}.$$

$$de=-\sum _{j=0}^n\frac{\partial g}{\partial q_j}d{q}_j+Td\eta -\frac{1}{{\rho }^2}\mathbf{\sigma }:{\mathbf{ϵ}}_{e}d\rho +\frac{1}{\rho }\mathbf{\sigma }:d{\mathbf{ϵ}}_e$$

$$de=\frac{1}{\rho }\mathbf{\sigma }:d{\mathbf{ϵ}}_e+Td\eta -\left(\frac{\partial g}{\partial \rho }+\frac{1}{{\rho }^2}\mathbf{\sigma }:{\mathbf{ϵ}}_e\right)d\rho -\sum _{j=1}^n\frac{\partial g}{\partial q_j}d{q}_j.$$

$$Q_0\equiv -\left(\rho \frac{\partial g}{\partial \rho }+\frac{1}{\rho }\mathbf{\sigma }:{\mathbf{ϵ}}_e\right)$$

$$Q_j\equiv -\rho \frac{\partial g}{\partial q_j}\text{(j=1, 2, 3, …)},$$

$$de=\frac{1}{\rho }\mathbf{\sigma }:d{\mathbf{ϵ}}_e+Td\eta +\frac{1}{\rho }\sum _{j=0}^n{Q}_{j}d{q}_j.$$

$$\frac{de}{dt}=T\frac{d\eta }{dt}+\frac{1}{\rho }\mathbf{\sigma }:{\dot{\mathbf{ϵ}}}_e+\frac{1}{\rho }\sum _{j=0}^n{Q}_j\frac{d{q}_j}{dt}.$$

$$\rho \frac{de}{dt}=\rho T\frac{d\eta }{dt}+\mathbf{\sigma }:{\dot{\mathbf{ϵ}}}_e+\sum _{j=0}^n{Q}_j\frac{d{q}_j}{dt}.$$

$$\rho \frac{de}{dt}=\mathbf{\sigma }:\nabla \mathbf{v}-\nabla \cdot \mathbf{q}+\rho s,$$

$$\mathbf{\sigma }:{\dot{\mathbf{ϵ}}}_e=\mathbf{\sigma }:\dot{\mathbf{ϵ}}-\mathbf{\sigma }:({\dot{\mathbf{ϵ}}}_p+{\dot{\mathbf{ϵ}}}_v).$$

$$\rho T\frac{d\eta }{dt}+\mathbf{\sigma }:\dot{\mathbf{ϵ}}-\mathbf{\sigma }:({\dot{\mathbf{ϵ}}}_p+{\dot{\mathbf{ϵ}}}_v)+\sum _{j=0}^n{Q}_j\frac{d{q}_j}{dt}-\mathbf{\sigma }:\dot{\mathbf{ϵ}}-\rho s+\nabla \cdot \mathbf{q}=0,$$

$$\rho T\frac{d\eta }{dt}=\mathbf{\sigma }:({\dot{\mathbf{ϵ}}}_p+{\dot{\mathbf{ϵ}}}_v)-\nabla \cdot \mathbf{q}+\rho s-\sum _{j=0}^n{Q}_j\frac{d{q}_j}{dt}.$$

$$d\eta =\frac{\partial \eta }{\partial \mathbf{\sigma }}:d\mathbf{\sigma }+\frac{\partial \eta }{\partial T}dT+\sum _{j=0}^n\frac{\partial \eta }{\partial q_j}d{q}_j.$$

$$d\eta =\frac{\partial \eta }{\partial \mathbf{\sigma }}:d\mathbf{\sigma }+\frac{\partial \eta }{\partial T}dT+\frac{\partial \eta }{\partial \rho }d\rho +\sum _{j=1}^n\frac{\partial \eta }{\partial q_j}d{q}_j.$$

$$\frac{d\eta }{dt}=\frac{\partial \eta }{\partial \mathbf{\sigma }}:\frac{d\mathbf{\sigma }}{dt}+\frac{\partial \eta }{\partial T}\frac{dT}{dt}+\frac{\partial \eta }{\partial \rho }\frac{d\rho }{dt}+\sum _{j=1}^n\frac{\partial \eta }{\partial q_j}\frac{d{q}_j}{dt}.$$

$$\frac{\partial \eta }{\partial \mathbf{\sigma }}=\frac{\partial }{\partial \mathbf{\sigma }}\frac{\partial g}{\partial T}=\frac{\partial }{\partial T}\frac{\partial g}{\partial \mathbf{\sigma }}=\frac{1}{\rho }\frac{\partial {\mathbf{ϵ}}_e}{\partial T}.$$

$$c_{\mathbf{\sigma }}\equiv \frac{\partial e}{\partial T}=\frac{\partial }{\partial T}\left(-g+T\eta +\frac{1}{\rho }\mathbf{\sigma }:{\mathbf{ϵ}}_e\right)=T\frac{\partial \eta }{\partial T}+\frac{1}{\rho }\mathbf{\sigma }:\frac{\partial {\mathbf{ϵ}}_e}{\partial T}.$$

$$\frac{\partial \eta }{\partial T}=\frac{1}{T}\left(c_{\mathbf{\sigma }}-\frac{1}{\rho }\mathbf{\sigma }:\frac{\partial {\mathbf{ϵ}}_e}{\partial T}\right).$$

$$\frac{\partial \eta }{\partial \rho }=\frac{\partial }{\partial \rho }\frac{\partial g}{\partial T}=\frac{\partial }{\partial T}\frac{\partial g}{\partial \rho }=-\frac{1}{\rho }\frac{\partial Q_0}{\partial T}-\frac{1}{{\rho }^2}\mathbf{\sigma }:\frac{\partial {\mathbf{ϵ}}_e}{\partial T}.$$

$$\frac{\partial \eta }{\partial q_j}=\frac{\partial }{\partial q_j}\frac{\partial g}{\partial T}=\frac{\partial }{\partial T}\frac{\partial g}{\partial q_j}=-\frac{1}{\rho }\frac{\partial Q_j}{\partial T},$$

$$\begin{array}{rl}\frac{d\eta }{dt}=& \frac{1}{\rho }\frac{\partial {\mathbf{ϵ}}_e}{\partial T}:\frac{d\mathbf{\sigma }}{dt}+\frac{1}{T}\left(c_{\mathbf{\sigma }}-\frac{1}{\rho }\mathbf{\sigma }:\frac{\partial {\mathbf{ϵ}}_e}{\partial T}\right)\frac{dT}{dt}\\ & +\left(-\frac{1}{\rho }\frac{\partial Q_0}{\partial T}-\frac{1}{{\rho }^2}\mathbf{\sigma }:\frac{\partial {\mathbf{ϵ}}_e}{\partial T}\right)\frac{d\rho }{dt}-\frac{1}{\rho }\sum _{j=1}^n\frac{\partial Q_j}{\partial T}\frac{d{q}_j}{dt}.\end{array}$$

$$\begin{array}{rl}\frac{d\eta }{dt}=\frac{1}{T}\left(c_{\mathbf{\sigma }}-\frac{1}{\rho }\mathbf{\sigma }:\frac{\partial {\mathbf{ϵ}}_e}{\partial T}\right)\frac{dT}{dt}+\frac{1}{\rho }\frac{\partial {\mathbf{ϵ}}_e}{\partial T}:\frac{d\mathbf{\sigma }}{dt}& -\frac{1}{{\rho }^2}\left(\mathbf{\sigma }:\frac{\partial {\mathbf{ϵ}}_e}{\partial T}\right)\frac{d\rho }{dt}\\ & -\frac{1}{\rho }\sum _{j=0}^n\frac{\partial Q_j}{\partial T}\frac{d{q}_j}{dt}.\end{array}$$

$$\begin{array}{r}\rho T\left[\frac{1}{T}\left(c_{\mathbf{\sigma }}-\frac{1}{\rho }\mathbf{\sigma }:\frac{\partial {\mathbf{ϵ}}_e}{\partial T}\right)\frac{dT}{dt}+\frac{1}{\rho }\frac{\partial {\mathbf{ϵ}}_e}{\partial T}:\left(\frac{d\mathbf{\sigma }}{dt}-\frac{1}{\rho }\mathbf{\sigma }\frac{d\rho }{dt}\right)-\frac{1}{\rho }\sum _{j=0}^n\frac{\partial Q_j}{\partial T}\frac{d{q}_j}{dt}\right]\\ -\mathbf{\sigma }:({\dot{\mathbf{ϵ}}}_p+{\dot{\mathbf{ϵ}}}_v)+\sum _{j=0}^n{Q}_j\frac{d{q}_j}{dt}-\rho s+\nabla \cdot \mathbf{q}=0,\end{array}$$

$$\begin{array}{rl}& \rho \left(c_{\mathbf{\sigma }}-\frac{1}{\rho }\mathbf{\sigma }:\frac{\partial {\mathbf{ϵ}}_e}{\partial T}\right)\frac{dT}{dt}+T\frac{\partial {\mathbf{ϵ}}_e}{\partial T}:\left(\frac{d\mathbf{\sigma }}{dt}-\frac{1}{\rho }\mathbf{\sigma }\frac{d\rho }{dt}\right)\\ & +\sum _{j=0}^n\left(Q_j-T\frac{\partial Q_j}{\partial T}\right)\frac{d{q}_j}{dt}-\mathbf{\sigma }:({\dot{\mathbf{ϵ}}}_p+{\dot{\mathbf{ϵ}}}_v)-\rho s+\nabla \cdot \mathbf{q}=0.\end{array}$$

$$\begin{array}{rl}\rho c_{\mathbf{\sigma }}\frac{dT}{dt}-\mathbf{\sigma }:{\alpha }_l\mathbf{I}\frac{dT}{dt}=& -\nabla \cdot \mathbf{q}+\rho s+\mathbf{\sigma }:({\dot{\mathbf{ϵ}}}_p+{\dot{\mathbf{ϵ}}}_v)\\ & -T{\alpha }_l\mathbf{I}:\left(\frac{d\mathbf{\sigma }}{dt}-\frac{1}{\rho }\frac{d\rho }{dt}\mathbf{\sigma }\right)\\ & -\sum _{j=0}^n\left(Q_j-T\frac{\partial Q_j}{\partial T}\right)\frac{d{q}_j}{dt}.\end{array}$$

$$\begin{array}{rl}(\rho c_p+p{\alpha }_v)\frac{dT}{dt}=& -\nabla \cdot \mathbf{q}+\rho s+\mathbf{\sigma }:({\dot{\mathbf{ϵ}}}_p+{\dot{\mathbf{ϵ}}}_v)+T{\alpha }_v\frac{dp}{dt}\\ & -pT\frac{{\alpha }_v}{\rho }\frac{d\rho }{dt}-\sum _{j=0}^n\left(Q_j-T\frac{\partial Q_j}{\partial T}\right)\frac{d{q}_j}{dt}.\end{array}$$

Full text

††journal: Tectonophysics

1 Introduction

The energy conservation principle can describe a wider range of geological and geodynamic phenomena when it takes into account energies involved with than the usual form involving only heat advection and diffusion because the general form can account for the energetics of deformation, which is often important in complicated phenomena. For instance, Hunt and Wadee [1991] show that non-periodic localized folding can be viewed as a superposition of folding modes that correspond to multiple local minima of the non-convex energy function and the non-convexity originates from deformational contribution to the system’s energy budget. On a similar note, Hobbs et al. [2011] propose that the feedback between shear heating and temperature-dependent viscosity can explain non-periodic and non-symmetric folding occurring in layers with small viscosity contrast, The classical Biot’s theory [e.g., Biot, 1961] predicts that folding would not occur in such a configuration. Influences of energy dissipation have been considered in the lithospheric scale as well. Regenauer-Lieb et al. [2006] consider the feedback between the energy dissipated due to inelastic deformation and changes in viscous and plastic material properties due to temperature changes resulting from heat converted from the dissipated energy. They find that the two-way feedback process can make the brittle-ductile transition (BDT) zone weaker than other parts of lithosphere although the classical strength envelopes predict the BDT zone of lithosphere to be the strongest [Ranalli and Murphy, 1987, Goetze and Evans, 1979, Brace and Kohlstedt, 1980]. Energy dissipated in the form of shear heating is shown to promote a necking in the subducting slab and ultimately lead to slab detachment [Gerya et al., 2004]. In the whole-mantle scale, Yuen et al. [1987] point out that feedback between rheology and dissipative energy in mantle convection is potentially an important mechanism that can warm up the mantle by several hundred degrees above the incompressible profile. Ita and King [1994] show that dissipative heating in the energy balance can significantly facilitate the vertical flow across the 660-km phase boundary.

However, computational long-term tectonic modeling, concerned about the long-term evolution of geological structures of various scales, has yet to fully embrace deformation energetics. Problems are centered around the fact that simplifications commonly made in this type of modeling preclude consistent thermo-mechanical coupling. For instance, energy balance is considered only in terms of heat advection and diffusion while thermo-mechanical feedback is realized only through temperature-dependent viscosity or shear heating.

More specifically, we identify three elements in kinematics and constitutive models that need to be incorporated into a thermodyanmic framework for long-term tectonic models. Firstly, we note that the elastic or plastic deformations are frequently assumed to be incompressible [e.g., Regenauer-Lieb and Yuen, 2003, Regenauer-Lieb et al., 2006, 2008, Connolly, 2009, Hobbs et al., 2011] although this assumption is neither required nor well-justified. Volumetric strain can have a significant effect on energy budget [Hunsche, 1991, Zinoviev and Ermakov, 1994] and is non-negligible during brittle deformations [e.g. Brace et al., 1966] and phase transformations Hyndman and Peacock [2003], Hetényi et al. [2011]. Secondly, thermal stresses are often ignored even though they can be a significant source of transient stresses and associated deformation [Choi et al., 2008, Choi and Gurnis, 2008, Korenaga, 2007]. A simple back-of-the-envelope calculation shows that a temperature change of 100 K in a perfectly confined rock body with a bulk modulus of 30 GPa can generate up to 90 MPa thermal stresses when the thermal expansion coefficient is $3\times 10^{-5}K^{-1}$. Although transient, thermal stresses of such a magnitude might be sufficient for driving permanent changes in state variables such as elastic damage or plastic strain under non-linear rheologies. Thirdly, strain weakening plasticity is sometimes considered as inconsistent with thermodynamic principles [e.g. Regenauer-Lieb et al., 2006]. However, frictional materials do show reduction 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 tectonic models concerned about brittle behaviors of rocks. In fact, the strain softening is perfectly legitimate in the light of the Clausius-Duhem inequality, a statement of the 2nd law of thermodynamics [e.g., Sec. 3.2 in Lubliner, 2008]. Rates and amounts of strain softening in rocks are still to be better constrained but reducing the uncertainty associated with them is an independent issue.

In this paper, we first derive a set of governing equations for thermo-mechanically coupled tectonic systems. We start from the generic thermodynamic principles closely following the procedure of  [Wright, 2002] to derive the governing equation that incorporate the three elements discussed above. We focus on the elasto-visco-plastic material that has compressible elasticity, thermal stress and strain wakening plasticity. Integration of such kinematic and constitutive models into the general energy balance will be crucial for realistically modeling geological systems. We then implement the derived governing equations in DES3D, an open source finite element code for geodynamic modeling, and verify the implementation semi-analytically. Finally, we explore the effects of thermo-mechanical coupling on the long-term evolution of large-offset normal faults with emphasis on the role of volumetric inelastic strain. Since extensibly studied and well understood, the normal fault systems would allow us to isolate the new effects introduced by the coupled physics.

2 Derivation of governing equations

2.1 Energy Balance Equation

Several theoretical works have derived a set of governing equations for a thermo-mechanical system from the general form of the thermodynamic principles. The common procedure is to relate the rate of change of the internal energy appearing in the statement of energy balance to that of thermodynamic potentials such as the Helmholtz free energy and the Gibbs free energy. Thermodynamic potentials involve 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 doing work at constant temperature” [p.263 in Malvern, 1969]. The main difference between the the Helmholtz and the Gibbs free energy is whether strain is an independent variable as in the former or stress is as in the latter. Since the definitions of these thermodynamic potentials involve the product of temperature and entropy, the energy balance principle takes an intermediate form involving the time derivative of entropy. The last step in deriving the temperature evolution equation is to express the time derivative of entropy in terms of that of temperature and other variables. Regenauer_Lieb_2003 start from the energy conservation principle stated in terms of the Helmholtz free energy to derive the partial differential equation for temperature evolution as well as other equations that are coupled with it (e.g., mass conservation equation and constitutive relations) for shear zone-developing systems in the Earth and other planets. Their final system of equations can consistently describe the feedback processes among energy, rheology and other variables such as grain size 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 energy conservation principle. Also starting from the energy balance equation in terms of the Helmholtz free energy, they include a non-dimensional variable quantifying the amount of damage along with temperature and infinitesimal elastic strain as independent variables of the free energy. By treating damage process as a source of entropy in the intermediate equation for entropy evolution, they derive a damage evolution equation that is proportional to the rate of free energy change with damage.

For completeness, we derive a temperature evolution equation from the generic energy balance principle involving the Gibbs free energy ($g$), following  Wright [2002]. This form of the energy balance principle is different from those of related studies in geodynamics  [Regenauer-Lieb and Yuen, 2003, Lyakhovsky et al., 1997] in that those earlier studies used another thermodynamic potential, the Helmholtz free energy. In the context of continuum thermodynamics, the Gibbs free energy is a function of Cauchy stress ($\boldsymbol{\sigma}$), temperature ($T$) and a set of internal variables ($q_{j}$, $j=0,1,2...n$) while the Helmholtz energy has elastic strain in place of Cauchy stress. The work by Wright [2002] inspired our choice of thermodynamic potential but one difference is that we assume infinitesimal strain while Wright [2002] used the finite strain kinematics.

Additive strain rate decomposition is allowed under our assumption of infinitesimal strain. Furthermore, the class of material we are interested in is elasto-visco-plastic. Under these assumptions, total strain rate is decomposed as follows:

[TABLE]

where $\dot{\boldsymbol{\epsilon}}_{e}$, $\dot{\boldsymbol{\epsilon}}_{p}$, and $\dot{\boldsymbol{\epsilon}}_{v}$ are elastic, plastic and viscous strain rate, respectively.

Following Wright [2002], we define the Gibbs free energy per unit mass as

[TABLE]

where $e$ is the internal energy per unit mass, $\eta$ is entropy per unit mass, $\rho$ is density and $\boldsymbol{\epsilon}_{e}$ is elastic strain tensor. Only the elastic strain appears in the definition because it is the only component of strain that can contribute to stored energy. Taking the totral differential of $g$ and rearranging terms, we get the differential of internal energy:

[TABLE]

Since $g$ = $g(\boldsymbol{\sigma},T,q_{j})$, the total differential of $g$ is also

[TABLE]

where $q_{0}=\rho$ according to the conventions of Wright [2002]. Elastic strain and entropy [c.f., Wright, 2002] are defined as

[TABLE]

From the equations (3), (4) and (5) we get

[TABLE]

Since $q_{0}=\rho$, the terms containing $dq_{0}$ and $d\rho$ can be grouped together. With the grouping, equation (6) becomes

[TABLE]

With two new notations

[TABLE]

and

[TABLE]

equation (6) is simplified to

[TABLE]

Differentiating the equation (9) with respect to time ($t$) we have

[TABLE]

Multiplying the equation (10) by $\rho$, we get the following equation for the time rate of change of internal energy per volume:

[TABLE]

To relate the material time derivative of the internal energy given in (11) to the energy balance principle, we recall the general form of the energy balance equation [e.g., Kennett and Bunge, 2008, Malvern, 1969, Wright, 2002]:

[TABLE]

where $s$ is a heat energy source or sink per mass, $\boldsymbol{q}$ is heat flux and $\boldsymbol{v}$ is velocity. According to the additive decomposition of strain rate in (1),

[TABLE]

Since the double dot product of the symmetric $\boldsymbol{\sigma}$ and the anti-symmetric part of $\nabla\boldsymbol{v}$ is zero, $\boldsymbol{\sigma}:\nabla\boldsymbol{v}=\boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}}$. By eliminating the time derivative of internal energy from equations (11) and (12) and then using (13), we get

[TABLE]

which is simplified to

[TABLE]

The final step to derive a partial differential equation for temperature from (14) is to relate the time derivative of entropy per mass to temperature. The equipresence principle  [Malvern, 1969] requires the entropy $\eta$ to have the same set of independent variables as the Gibbs free energy. In other words, the entropy per unit mass is also a function of the Cauchy stress ($\boldsymbol{\sigma}$), the temperature ($T$) and a set of internal variables ($q_{j}$, $j=0,1,2...n$). The total differential of $\eta=\eta(\boldsymbol{\sigma},T,q_{j})$ is

[TABLE]

Identifying the first internal variable ($q_{0}$) with density ($\rho$) again, we get

[TABLE]

When differentiated with respect to time, the equation (15) becomes

[TABLE]

According to (5),

[TABLE]

From the definition of the specific heat at constant stress, $c_{\boldsymbol{\sigma}}$, we get the following identity:

[TABLE]

Using (18), we get

[TABLE]

For convenience, we express partial derivatives of entropy per mass with respect to density and other internal variables in terms of $Q_{0}$ and $Q_{j}$ ($j=1,2,3,\ldots$), which are defined in (7) and (8).

[TABLE]

and

[TABLE]

for $j=1,2,3,\ldots$. Plugging (17), (19), (20) and (21) and into the equation (16), we get

[TABLE]

This equation can be further simplified to

[TABLE]

Substituting the equation (22) into the equation (14) we get:

[TABLE]

which is simplified to

[TABLE]

Strain due to temperature change is assumed to be isotropic: i.e., $\partial\boldsymbol{\epsilon}_{e}/\partial T=\alpha_{l}\boldsymbol{I}$, where $\boldsymbol{I}$ is the identity tensor and $\alpha_{l}$ is linear thermal expansion coefficient and 1/3 of the volumetric thermal expansion coefficient, $\alpha_{v}$. Under this assumption, the equation (23) is rearranged as follows:

[TABLE]

Using the identity $\boldsymbol{\sigma}:\alpha_{l}I=-3\,\alpha_{l}\,p=-\alpha_{v}\,p$ where $p=-\frac{1}{3}\sigma_{ii}$, we have the following final form of the energy balance equation:

[TABLE]

The last term of (25) corresponds to the changes in energy due to internal variables only. These terms are often parametrized into a coefficient for the inelastic power term as in $\chi\boldsymbol{\sigma}:(\dot{\boldsymbol{\epsilon}}_{p}+\dot{\boldsymbol{\epsilon}}_{v})$  [Regenauer-Lieb and Yuen, 2003, Wright, 2002]. With this parametrization, the equation (25) is simplied to

[TABLE]

$\chi$ close to 1 means that the energy change due to changes in internal variables is negligible relative to inelastic power and most of the inelastic power is converted to heat production. To our knowledge, the value of $\chi$ is not very well constrained for rocks but at least for metals, it is either close to 1 or saturates towards 1 with increasing plastic strain [Wright, 2002]. Regenauer-Lieb and Yuen [2003] used 0.85 as the value of $\chi$. $\chi$ is assumed to be 1 in this study.

2.2 Mass balance equation

The form of the mass conservation equation involving the material derivative is

[TABLE]

In the Lagrangian description of motion that we are going to adopt, the above time derivative is understood as the partial derivative with respect to time, not as the material time derivative. We substitute (27) into the last term of (26) to get

[TABLE]

2.3 Thermoelastic constitutive equations

Assuming the linear isotropic elasticity and temperature change $\delta T$, we use the following thermoelastic constitutive equations [e.g., Boley and Weiner, 1997]:

[TABLE]

where $\lambda$ and $G$ are the Lamé’s constant and $K$ is the bulk modulus defined as $\lambda+2G/3$. These thermal elastic constitutive equations become the basis for viscoelastic or elastoplastic constitutive models as described in [Choi et al., 2013b].

2.4 On the use of non-associated strain-softening plasticity

Using a non-associated, strain-weakening plasticity for inducing shear bands as in [Choi et al., 2013b] has been described as if it violated the 2nd law of thermodynamics [e.g., Regenauer-Lieb and Yuen, 2003, Regenauer-Lieb et al., 2006, 2008]. An accompanying criticism on the strain-weakening plasticity is that the softening rules often used in tectonic modeling lack reliable constraints. While agreeing that the lack of experimental and observational constraints on the strain-softening rules is problematic, we would like to point out that this problem does not invalidate the approach itself. In fact, shear band formation induced by strain softening is perfectly legitimate in the light of the 2nd law of thermodynamics. The source of previous confusion must be in the erroneous notion that only a strain hardening plasticity is “stable” in the sense that the plastic dissipation is positive therefore the 2nd law of thermodyanmics is satisfied. However, neither strain softening nor non-associated flow rule necessarily violates the maximum plastic dissipation postulate as shown by  Lubliner [2008, see Sec. 3.2.]

The maximum plastic dissipation postulate requires that the yield surface be convex in the principal stress space and the isotropic strain-softening plasticity model often used in tectonic modeling! satisfy this requirement by maintaining the convexity. For instance, when a strain-softening in the Mohr-Coulomb plastic model is realized through reduction in cohesion and friction angle, the hexagonal cone-shaped, therefore convex, yield surface in the principal stress space always remain convex except in the degenerate case where both friction angle and cohesion are zero.

Under a typical setting for numerical tectonic modeling in which plastic flow is incompressible, i.e., dilation angle is 0*∘, and the friction angle is about $30^{\circ}$, the principal 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 mean that such a non-associated flow rule violates the maximum plastic dissipation. Friction angles usually assumed in tectonic models are never much greater than 30∘* and dilation angles are less than that. Since principal stress and plastic strain orientations are normal to a yield surface and a flow potential, respectively, the difference between 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 definition of plastic dissipation rate is equivalent to the projection operation, an acute angle between the principal stresses and plastic strain rates guarantees a positive plastic dissipation that indicates the satisfaction of the maximum plastic dissipation postulate.

These considerations validate our adoption of the strain-weakening, non-associated, Mohr-Coulomb plastic model.

3 Benchmarks

We verify the governing equations implemented into DES3D, an unstructured finite element code for long-term tectonic modeling [Choi et al., 2013b], using three benchmark problems. The benchmarks are derived from the standard oedometer test [e.g., Davis and Selvadurai, 2002, Choi et al., 2013b]. A cubic block of the Mohr-Coulomb plastic material is compressed in one direction while motions in the other directions are restricted (Fig. 1). The symmetry and boundary conditions of the problem make it sufficient to discretize a square domain with two triangular elements. The simplicity of the problem allows for at least semi-analytic solutions. In benchmark-1, we include only the plastic power $(\boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}}_{p})$ as a heat source and ignore the diffusion term. The density is assumed to be constant. Benchmark-2 solves the following equation

[TABLE]

and the density is updated according to eq. (27). Benchmark-3 verifies the thermal stress calculations in DES3D for a uniform temperature that increases linearly in time. We intentionally use a low density (1 kg/m3) to make the contribution from plastic power non-negligible. Parameters used in the benchmarks are listed in Table 1.

3.1 Analytic solutions for the oedometer test

3.1.1 Total strain rate and strain

The model geometry and boundary conditions give the following components of total strain rate ($\dot{\boldsymbol{\varepsilon}}$):

[TABLE]

where $t$ is time, $v_{x}$ is the boundary velocity set to be $-10^{-5}$ m/s and $L$ is the edge length of the cube equal to 1 m.

The components of total strain ($\boldsymbol{\varepsilon}$) are given as

[TABLE]

3.1.2 Pre-yielding stress

Before yielding, stresses are updated according to the linear isotropic elasticity. In terms of the Lamé’s constants ($\lambda$, $G$), stress components are given as:

[TABLE]

3.1.3 Mohr-Coulomb yield function and flow potential

We use the following forms of the Mohr-Coulomb yield function ($f$) and flow potential ($g$):

[TABLE]

where $N_{\phi}$ is defined as $(1+\sin\phi)/(1-\sin\phi)$, $\phi$ is the firction angle, $\psi$ is the dilation angle and $C$ is the cohesion.

3.1.4 Yielding time ($t_{cr}$)

We denote the time when yielding occurs for the first time as $t_{cr}$. Since $f(\sigma_{xx},\sigma_{yy})=0$ at $t=t_{cr}$,

[TABLE]

Solving the above equation for $t_{cr}$, we get

[TABLE]

3.1.5 Plastic strain rate and strain

We get the following expressions for plastic strain for $t\geq t_{cr}$:

[TABLE]

where $\beta(t)$ is the plastic consistency paramter to be determined. As in the classical plasticity theory [e.g., Lubliner, 2008], the consistency parameter is zero before yielding ($t\leq t_{cr}$).

The plastic strain rates are similarly defined in terms of the rate of consistency parameter ($\dot{\beta}(t)$):

[TABLE]

3.1.6 Elastic strain after yielding

The above definitions of total and plastic strain lead to the following expressions for elastic strain as a function of time:

[TABLE]

3.1.7 Post-yielding stress and determination of $\beta(t)$

After yielding occurs, stresses are updated as follows:

[TABLE]

Plugging (49) to (51) into the above equations, we get

[TABLE]

For $t\geq t_{cr}$, these stresses should always satisfy the yield condition:

[TABLE]

Substituting (56) and (57) for the stress components in the yield condition, we get

[TABLE]

Solving the above equation for $\beta(t)$, we get

[TABLE]

The time derivative of $\beta(t)$ is given as

[TABLE]

3.1.8 Pressure and pressure rate

Since $p=-(\sigma_{xx}+\sigma_{yy}+\sigma_{zz})/3$, we get the following from (56), (57) and (58):

[TABLE]

where $K$ is the bulk modulus defined as $\lambda+2/3G$.

3.1.9 Volume change rate

[TABLE]

3.1.10 Density

Using the expression for $\nabla\cdot\mathbf{v}$, the mass balance equation becomes

[TABLE]

Solving for $\rho(t)$, we get

[TABLE]

where $\rho_{0}$ is the reference density at $t=0$.

3.2 Benchmark-1 : Plastic power only, no thermal stress

The equation we are going to solve is

[TABLE]

where

[TABLE]

Since $\rho$ and $c_{p}$ are assumed to be constant,

[TABLE]

We use the forward Euler scheme to integrate the above time derivative of $T$ with the initial condition $T(0)=273$ K.

Results of the test are shown in Fig. 2. The relative error of the DES3D solution relative to the semi-analytic one is 0.004 %.

3.2.1 Benchmark-2 : The full equation, no thermal stress

The equation to solve is

[TABLE]

Using the derived expressions for the involved quantites and the forward Euler scheme, we can integrate the above equation for the initial condition, $T(0)=273.0$ K.

Results of the test are shown in Fig. 3. The relative errors of the temperature and density from DES3D are 0.01 % and 0.00003 %, respectively.

3.3 Benchmark-3: Thermal stresses under a prescribed temperature change

We verify the implementation of the thermoelastic constitutive equations in DES3D. For simplicity, we prescribe a uniform temperature field, which is initially 273 K and increases linearly in time as

[TABLE]

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 solutions for the oedometer test involving thermal stress. In the elastic regime, the stresses are given as

[TABLE]

Using the expressions, (56), (57) and (58), we get

[TABLE]

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 in terms of elementary functions. Following the procedure for determining $\beta(t)$ and $\dot{\beta}(t)$ described in Sec. 3.1.7, we get

[TABLE]

and

[TABLE]

As before, $\beta(t)=\dot{\beta}(t)=0$ for $t<t_{cr}$.

Differential stress arising in this benchmark is the same with those of the isothermal case but pressure in this test is greater due to the compressional thermal stress caused by the prescribed temperature increase. As a result, the first yielding in benchmark-3 should occur at a greater value of differential stress, or equivalently, displacement than in the isothermal case. In other words, than in the isothermal case. The results of benchmark-3 shows in Fig. 4 are consistent with this expectation. The relative errors of $\sigma_{xx}$ and $\sigma_{yy}$ computed with DES3D are 0.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 the thermo-mechanically coupled governing equations on a more geologically relevant problem. In this study, we choose the large-offset normal fault evolution [e.g. Lavier et al., 2000] as an example. Studied extensively and understood well, this system will facilitate identification and attribution of differences between the models of the newly-introduced physics and the uncoupled, isothermal ones.

We create a normal fault in an extending Mohr-Coulomb elasto-plastic layer which is initially 100 km long and 10 km thick (Fig. 5). Both sides of the layer are pulled at a constant velocity of 0.5 cm/yr. The bottom boundary is supported by the Winkler foundation [Watts, 2001, pp.95]. To induce strain localization, we decrease 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 from there. We impose topographic smoothing of the diffusion type with a transport coefficient of $10^{-7}m^{2}/s$ [Turcotte and Schubert, 2014, pp. 225].

We set up five models that differ in the form of the energy balance equation and in the presence of volumetric plastic strain. In model-1, we do not consider any heat-generating mechanism. As a result, the initial temperature, which is uniformly $0^{\circ}$C, does not change in time. The behavior of model-1 should be similar to that of the “unlimited, footwall-snapping” mode of [Lavier et al., 2000]. In the rest of models (model-2 to 5) we solve the full energy balance equation (28) with the heat source/sink term, $s$, equal to zero and with the temperature feedback to rheology through thermo-elasticity. To explore the effects of volumetric plastic strain, we use a non-zero dilation angle in model-3 to 5. The dilation angle is fixed in model-3 while we reduce it to zero as the accumulated plastic strain increases to a prescribed value of 1.0 in model-4 and 5. The model-5 is the same as model-4 but we intentionally include only the deviatoric part of the plastic power ($\boldsymbol{\sigma}:\dot{\boldsymbol{\varepsilon}}_{p}$) in equation (28). By comparing model-4 and 5, we can decide whether the volumetric plastic power can be ignored. Table 2 summarizes the differences among the models. Table 3 shows the list of the parameters used in the models.

We see noticeable differences in fault geometry and shape of the core complex between model-1 and model-2 that solves the full energy balance equation and includes thermal stresses. The overall behavior of the faults in model-1 is similar to the unlimited, footwall snapping mode of Lavier et al. [2000]. The geometry of the primary fault of model-2 is almost the same as that of the model-1 until about 17 km of extension. Around this time, however, model-2 forms a secondary fault but model-1 does not yet. After 20 km of extension, more secondary faults start form in both models but their geometry and location are not identical. For instance, a fault block forms in model-2 after about 30 km of extension when a new high-angle fault forms next to the first secondary fault. Model-1 does not develop a corresponding structure including a fault block after the same amount of extension. Model-1 forms a single secondary fault initially, which eventually connects with the second secondary fault after about 40 km of extension. At 30 km of extension model-2 has three well-developed secondary faults while model-1 has two. Initially, the shape of core complex remains same until 10 km of extension. Because of different fault geometry the shape of the core complex started to differ from each other from 20 km of extension. At 40 km of extension the core-complex of the model-2 becomes elongated with the almost flat surface while the model-1 has more rounded one.

Figure 7 shows the distribution of temperature in model-2 and the corresponding thermal stresses after 10 and 30 km of extension. The temperature increase near the active faults are greater than 15 *∘*C and the magnitude of thermal stress is as large as 40 MPa. The thermal stress magnitude is a non-negligible fraction, about 10 %, of lithostatic stress near the bottom. Thus, the different behaviors of the models can be attributed to the presence of thermal stresses only in model-2.

Identical with model-2 except a non-zero dilation angle of 10*∘, model-3 develops an unrealistically elevated core-complex with a relief greater than 10 km and shear bands that are much wider than those in model-1 and 2. Expansion of shear bands when dilation angle is non-zero is kinematically expected. The dilation angle fixed at 10∘* sustains a higher level of plastic power than the non-dilational cases zero dilation angle. The greater plastic power leads to a greater amount of temperature change and thermal stresses. The thermal stresses pushes up the free-traction top surface, creating the highly elevated core-complex.

To isolate the effect of volumetric plastic power, we construct the model-4 and 5 that are identical except that model-4 considers the total plastic power while model-5 includes only the deviatoric part of the plastic power. In order to prevent the unrealistic expansion of shear bands as in model-3, the dilation angle is reduced with accumulated plastic strain in both model-4 and 5 [Choi and Petersen, 2015].

Models-4 and 5 do not show notable differences until 10 km extension (Figure 9). However, after about 20 km of extension, they start to exhibit differences in the fault geometry as well as in the timing of rider block formation. A rider block starts to form at 17 km of extension in model-4 but model-5 forms one only after more than 37 km of extension. After 40 km of extension, model-5 has three secondary faults produced by the footwall snapping but model-4 still has only two secondary fault. The topographic relief in model-4 and 5 is about 5 km, which is much greater than the relief of the previous models for the large offset normal fault [Lavier et al., 2000, Choi et al., 2013a] but about half of that of model-3 with a fixed dilation angle 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 greater than that of this study to have 1-2 km topographic relief as in the previous studies on large-offset normal faults. The noticeable differences between model-4 and 5 suggest to include the total plastic power in the energy balance rather than only the deviatoric component. Ignoring the volumetric plastic power as in model-5 is inconsistent with the kinematics in the first place.

5 Conclusions

We derive a temperature evolution equation based on the energy balance principle that accounts for deformation-related energy changes as well as the conventional heat energy diffusion and advection. An explicit-time finite element solution procedure for the derived equation is implemented in DES3D, an unstructured mesh finite element solver for geodynamics. We verify the implementation using three benchmark tests that have semi-analytic solutions. Benchmark-1 includes only the total plastic power term in the energy balance equation. Benchmark-2 includes the full energy equation as well as the mass balance equation. Benchmark-3 prescribes a temperature change linear in time and computes thermal stresses. The numerical solutions from DES3D for all the benchmarks show an excellent match with the corresponding semi-analytic solutions. Coupling the verified implementation of the full temperature evolution with the strain weakening Mohr-Coulomb plastic rheology through thermal stresses, we also explore the long-term behaviors of a large-offset normal fault. We find that the temperature arising mostly from the inelastic power term in the full energy balance has non-negligible effects on the long-term evolution of normal fault. For instance, when the plastic power is considered, temperature increases by more than 15*∘*C and the magnitude of the associated thermal stress is as large as 40 MPa. These extra stresses appear to promote the formation of secondary faults and an elongated core complex. When the dilational plastic deformation is enabled, our models develop even more differences in faulting behaviors such as the formation of a rider block and the great topographic relief, compared to the previous models for the large-offset normal fault with the same parameters and geometry. We conclude that the new coupled governing equations presented in this study has significant impact on long-term tectonics. Although not so strong as to cause a fundamental shift in faulting modes under tested conditions, the effects of deformation energetics might be critically important in other geologically relevant systems and conditions.