THMC instabilities and the Multiphysics of Earthquakes

TL;DR

This paper explores how coupling Thermo-Hydro-Mechano-Chemical processes can lead to volumetric instabilities that potentially trigger earthquake shear instabilities through resonance phenomena in a reaction-diffusion framework.

Contribution

It introduces a novel approach analyzing THMC coupling and cross-diffusivities as a mechanism for earthquake precursor instabilities.

Findings

Cross-diffusivities induce slow pressure waves.

Resonance in long-wavelength limit triggers seismic events.

Volumetric instabilities may precede shear instabilities.

Abstract

Theoretical approaches to earthquake instabilities propose shear dominated instabilities as a source mechanism. Here we take a fresh look at the role of possible volumetric instabilities preceding a shear instability. We investigate the phenomena that may prepare earthquake instabilities using coupling of Thermo-Hydro-Mechano-Chemical (THMC) reaction-diffusion equations in a THMC diffusion matrix. We show that the off-diagonal cross-diffusivities can give rise to cross-diffusion slow (no inertia) pressure waves and propose that a resonance phenomenon of the long-wavelength limit causes constructive interference and triggers the seismic moment release.

Equations

$$\frac{\mathcal{D}{C}_i}{\mathcal{D}t}=D_{ii}{\nabla }^2{C}_i+\sum _{j\ne i}{D}_{ij}{\nabla }^2{C}_j+r_i.$$

$$\frac{\mathcal{D}{C}_i}{\mathcal{D}t}=\sum _{k=1}^N{D}_{ij}{\nabla }^2{C}_k+r_i,$$

$${{\mathbf{D}}_{ij}}=\begin{bmatrix}{{\color[rgb]{1,0,0}{D}_{TT}}}&{{L}_{TH}}&{\color[rgb]{0,0,1}{{L}_{TM}}}&{{L}_{TC}}\\ {{L}_{HT}}&{\color[rgb]{1,0,0}{{D}_{HH}}}&{{L}_{HM}}&{{L}_{HC}}\\ {\color[rgb]{0,0,1}{{L}_{MT}}}&{{L}_{MH}}&{\color[rgb]{1,0,0}{{D}_{MM}}}&{{L}_{MC}}\\ {{L}_{CT}}&{{L}_{CH}}&{{L}_{CM}}&{\color[rgb]{1,0,0}{{D}_{CC}}}\\ \end{bmatrix}.$$

$$\frac{\mathcal{D}{C}_1}{\mathcal{D}t}=D_{HH}\frac{{\partial }^2{C}_1}{\partial x^2}+D_{HM}\frac{{\partial }^2{C}_2}{\partial x^2}$$

$$\frac{\mathcal{D}{C}_2}{\mathcal{D}t}=D_{MM}\frac{{\partial }^2{C}_2}{\partial x^2}+D_{MH}\frac{{\partial }^2{C}_1}{\partial x^2}-\lambda C_2$$

$$\left[\begin{array}{c}{C}_1\\ C_2\end{array}\right]=\left[\begin{array}{c}{C}_1\exp \left(ikx+\omega t\right)\\ C_2\exp \left(ikx+\omega t\right)\end{array}\right],$$

$$(D_{HH}-D_{MM}-\frac{D_{MM}}{D_{HH}{k}^2}{)}^2+4D_{HM}{D}_{MH}<0.$$

$$(D_{HH}-D_{MM}-1{)}^2+4D_{HM}{D}_{MH}<0.$$

Full text

I Introduction

We present the hypothesis that the universality of cross-diffusion phenomena based on multiscale, multiphysics coupling underpins the scale invariance of the Gutenberg, Richter frequency-magnitude relationship of earthquakes (Sethna et al., 2001). We propose a new mathematical approach for coupling the rates of chemical reactions, failure and fluid flow instabilities and thermo-mechanical instabilities of materials. The approach gives theoretical insights into the processes that are commonly described by empirical relationships.

In our multiphysics approach to earthquakes we consider local chemical reactions and mass exchange processes as the smallest thermodynamic driving force triggering a thermodynamic flux. The approach emphasises the important role of emergent patterns (microstructure) triggering energy-geometry interactions that cascade through space and time via the link of THMC cross-diffusion coefficients introduced in this paper. These coefficients can be activated in far from equilibrium conditions, which allows coupling between the scales. The important factor of cross-scale coupling by cross-diffusion was overlooked in earlier reviews of multiscale coupling and multiphysics approaches (Regenauer-Lieb et al., 2013a, b) to earth physics problems. We propose here that earthquakes can appear if the large-scale geodynamic driver fires the full cascade of THMC cross-diffusion coefficients through reaching a critical magnitude of reactive THMC source terms.

II Multiscale cross-diffusion model

The conventional definition of thermodynamic forces and fluxes in a THMC-system is given in Table 1. Thermodynamic forces are understood as the spatial gradients of a THMC-system. They drive time gradient thermodynamic fluxes. At the molecular scale, the thermodynamic force is, for instance, the gradient of a chemical species, which drives a chemical flux (C) parameterized by Fick’s law, where ${{D}_{C}}$ is the chemical diffusivity. Such chemical reactions can in turn become a driver for the next scale up in the multiscale system. For example, a chemical fluid release reaction in an HMC-coupled system can create mechanical pressure gradients due to the associated volume change and fluid release/precipitation, which in turn induces mechanical (M) stress diffusion ${{D}_{M}}$ with a Beltrami-Michell-type compatibility condition (Rice and Cleary, 1976) and Darcy-type fluid (H) diffusion ${{D}_{H}}$. In the mechanical sense, the chemical reaction creates an overpressure ${{\bar{p}}_{M}}$ in the matrix (Hu et al., 2019). The endo-/exo-thermal chemical reaction also induces a gradient in temperature which in turn drives a thermal Fourier flux (T), with thermal diffusivity ${{D}_{T}}$.

The afore mentioned cross-scale coupling introduces source/sink terms in the individual conservation laws identified by ${{r}_{T}},{{r}_{H}},{{r}_{M}}$ and ${{r}_{C}}$ respectively. The conservation laws must also be extended to allow a multitude of THMC processes to occur simultaneously, which can introduce cross-diffusion fluxes.

II.1 THMC cross-diffusion model formulation

In a chemical system with two species, for instance, cross-diffusion is the phenomenon, in which a flux of species $A$ is induced by a gradient of species $B$ (Vanag and Epstein, 2009). In more general THMC terms, cross-diffusion is the phenomenon where a gradient of one generalised thermodynamic force of species ${C}_{j}$ drives another generalised thermodynamic flux of species ${C}_{i}$, described by (Manning, 1970)

[TABLE]

The species $j$ is identified as the cross diffusion species other than the species $i$ . Introducing a fully populated ($N\times N$) diffusion matrix ${{\mathbf{D}}_{ij}}$, equation (1) can also be written as

[TABLE]

The diffusion matrix ${{\mathbf{D}}_{ij}}$ for the equivalent THMC processes is

[TABLE]

The classical (self-)diffusive length scale of each THMC process is defined by $\sqrt{D_{ii}t}$. The reactive source source term $r_{i}$, corresponding to each thermodynamic flux ${C}_{i}$, can be zero or non-zero. The approach presented here introduces a new cross-scale coupling length scale that provides a link between the self-diffusion length scales of the THMC processes. Although not essential to the formulation we propose that this cross-diffusive length scale $\sqrt{D_{ij}t}$ is quantized and defined by discrete material properties which can be evaluated by measuring the velocity of cross-diffusion waves (Regenauer-Lieb et al., 2019). These (quantized) cross-diffusion coefficients link the gradient of a thermodynamic force ${C}_{j}$ of one THMC process to the flux of another kind.

Under such circumstances, the criterion for nucleation of diffusion waves is fully characterized by the diffusion matrix ${{\mathbf{D}}_{ij}}$ with the diagonal elements ${{D}_{ii}}$ describing the normal self-diffusion and the off-diagonal elements the cross-diffusion processes enabling coupling across scales. For consistency with the second law, all eigenvalues of the diffusion matrix must be real and positive, and hence the determinant of the diffusion matrix $\det\left(\mathbf{D}_{ij}\right)>0$ as well as the trace of the diffusion matrix must be larger than zero.

The criterion for instability can now be expressed by two matrices, the self- and cross-diffusion matrices, respectively. The cross-diffusion matrix must satisfy the Onsager-Casimir symmetry condition which acknowledges antisymmetry ${{D}_{ij}}=\pm{{D}_{ji}}$ (Onsager, 1931; Casimir, 1945). Onsager (Onsager, 1931) assumed that the parameters characterizing the system have identical sign in their kinetic coefficients while Casimir (Casimir, 1945) extended Onsager’s considerations to include kinetic parameters that have opposite sign as can be the case for multiphysics cross-diffusion coupling. Any local perturbation that leads to complex eigenvalues of ${\mathbf{D}_{ij}}$ must lead to oscillatory instabilities relaxing to the equilibrium state (Vanag and Epstein, 2009). This statement is at the heart of the nucleation of the cross-diffusional wave phenomenon. We will show in the following sections that the above described cross-diffusion formulation (Vanag and Epstein, 2009) can be extended to develop a generic multiphysics and multiscale THMC coupling approach to earthquake instabilities. This formulation implies a coupled cascade of THMC feedbacks over multiple diffusional length scales honouring the reciprocal multiscale interplay of thermodynamic forces and fluxes.

II.2 Cross-diffusion and its role in coupling of instabilities across scales

By extending the diagonal diffusion matrix through the cross-diffusion coefficients in equation (2), a new class of mechanical instabilities is revealed that incorporates chemical instabilities at its lowest scale (Regenauer-Lieb et al., 2019). For the simple case of hydro-mechanical (HM) coupling we have recently reported discovery of the cross-diffusional pressure wave phenomenon (Hu et al., 2019).

The important element of cross-diffusion waves for earthquake physics is their capability to link one thermodynamic force with a thermodynamic flux at a different scale, thus synchronising the dynamics over vastly different diffusional time and length scales. This important aspect of earthquake physics was previously overlooked. The approach raises the possibility that dissipative waves can be detected prior to earthquake instabilities.

III Application to earthquake physics

Theoretical approaches to the physics of earthquake instabilities originally were conceived as a shear instability on a frictional sliding surface (Brace and Byerlee, 1966). The role of pressure on the dynamics of the slider was derived empirically by laboratory experiments defining a rate-and-state variable friction law (Dieterich, 1972, 1987; Ruina, 1983; Tse and Rice, 1986; Rice et al., 2001). Instabilities through shear heating feedback were later considered to play an important role (Ogawa, 1987; Regenauer-Lieb and Yuen, 1998; Braeck and Podladchikov, 2007). Additionally, a thermally induced fluid pressurization term was found to be an important component for accelerated creep (Vardoulakis, 2001; Rice, 2006). Another important ingredient of the earthquake instability was thought to be the coupling of scales, where at least two different processes, operating at different time and length scales, interact (Ohnaka, 2003). The approach presented here summarizes these effects in the matrix ${{\mathbf{D}}_{ij}}$ and enables upscaling via renormalization group theory (Regenauer-Lieb et al., 2013a; Hanasoge et al., 2017).

III.1 Cross-diffusion waves

A fully populated diffusion matrix provides the opportunity to extend the postulate of coupling different scale processes in earthquake mechanics (Ohnaka, 2003). We present a high-level application of the basic principles in order to transfer knowledge from other disciplines to this new field of research (e.g. mathematical biology, computational chemistry, quantum optics). We use closed form solutions developed in these fields, supplemented by our own solutions (Regenauer-Lieb et al., 2013b; Veveakis and Regenauer-Lieb, 2015), for schematic illustrations and in-depth discussions of how the cross-diffusion concept may be applied to the earthquake phenomenon.

Exact analytical propagating wave solutions of the cross-diffusion wave phenomenon are rare. The FitzHugh-Nagumo oscillator (Zemskov et al., 2017) is by far the best analysed prototype of an excitable system with short lived-spikes and a slow recovery. The FitzHugh-Nagumo cross-diffusion model has originally been developed to describe the propagation of electrical impulses in nerves. We use a generic application of this oscillator to illustrate the potential feedback of just two THMC processes and the conditions that prepare the formation of a future fault plane. The real waveforms of the transient cross-diffusion waves are expected to be more complex.

III.1.1 Cross-diffusion wavenumbers

The wavenumbers of the cross-diffusion waves are controlled by the competition of the diffusion processes determining the diffusion length and number of possible oscillations per unit length. Of particular interest for the earthquake instabilities is the resonant mode with the longest possible wavelength.

Below, we discuss the results of a 1-D linear stability analysis for just two coupled reaction-diffusion equations to illustrate the role of the wavelength (or its inverse the wavenumber) for the nucleation of cross-diffusion waves. Laboratory evidence of cross-diffusion waves with HM coupling and a detailed derivation of the equations can be found in ref. Hu et al. (2019). For HM coupling, we set thermal and chemical self- and cross-diffusion coefficients to zero in equation (2) and using a single source term characterized by the ratio of the self-diffusivities of the two processes:

[TABLE]

where $\lambda={D_{MM}}/{D_{HH}}$, is a non-dimensional parameter denoting the ratio of mechanical diffusivity over the fluid diffusivity. The response of the system to wave perturbations is tested by applying a small plane wave perturbation of the type

[TABLE]

where $k$ is the wavenumber and $\omega$ the angular frequency. The perturbation yields condition for instability when (Hu et al., 2019)

[TABLE]

A cross-diffusional wave can nucleate if ${{D}_{HM}}$ and ${{D}_{MH}}$ are nonzero, of opposite sign and sufficiently large to overcome the self-diffusion processes in the first term, which is always positive. The instability criterion also suggests a critical lower limit $k_{min}$ (longest wavelength) for which no cross-diffusion wave can nucleate as the first term becomes dominant. This minimum wavenumber (longest wavelength) can be constrained further by the long wavelength limit of the system size. We argue that this size is a fraction of the system size, defined by long-wavelength Korteweg-deVries equation solution (Regenauer-Lieb et al., 2013b; Veveakis and Regenauer-Lieb, 2015) where the cross-diffusion coefficients are adiabatically eliminated. This long-wave soliton solution allows constructive interference of the cross-diffusional waves in normal resonance mode.

III.1.2 Resonance on the lowest possible wavenumber

In the context of an earthquake instability, shorter wavelength cross-diffusional waves might lead to a small dispersion limit of a coordinated long-wavelength instability. This resonance mode may arise through incident cross-diffusional waves colliding with reflected cross-diffusional waves in a finite width around the future fault core. On the grounds of the above theoretical considerations we postulate a series of evolutionary steps that may lead to this macro-scale instability. These are:

• •

Step 1 (illustrated in Figure 1); Upon ongoing geodynamic loading a zone may form, where the rock is loaded past the elastic limit, and the plastic yield stress of the rock is reached. This defines the system size of reflecting cross-diffusion waves.

• •

Step 2 (illustrated in Figure 2); If the criterion in equation (6) is met, cross-diffusional waves are expected to initially bounce back and forth inside the plastically deforming domain as wave packets with dominant finite wavenumbers (Biktashev and Tsyganov, 2016; Tsyganov and Biktashev, 2014), described in the next chapter on waveforms. In the hydromechanically coupled problem (equation 4), these waves are propagating fronts of internal volumetric fluid-solid mass exchange and appear as local fluid overpressure. Figure (2) shows conditions that may lead to the formation of the future fault plane by collision of cross-diffusion waves approaching each other from opposite sides of the fault plane. Collisions of cross-diffusion waves can lead to accumulation of material damage discussed in the section on waveforms.

• •

Step 3 (illustrated in Figure 3); The long wavelength limit of cross-diffusional waves is defined by a fraction of the width of the plastic zone. A large-scale instability is then understood as the longest wavelength normal (resonance) mode to form a standing wave. Earthquakes may occur prior to this situation when the volumetric damage accumulated in Step 2 lubricates a shear instability. For extremely slow deformation the purely volumetric long-wavelength solution of the standing P-wave can fully develop and is shown in Figure 3. We use the analytical solution of the Korteweg-de Vries equation (Regenauer-Lieb et al., 2013b; Veveakis and Regenauer-Lieb, 2015) and emphasize that the solutions for Figures 2 and 3 are for illustrative purposes only.

The earthquake trigger for positive interference necessary for the appearance of a normal mode standing wave is defined by substituting $k_{min}=\sqrt{\frac{D_{MM}}{D_{HH}}}$, obtained from the standing wave-limit, into the instability criterion (6):

[TABLE]

This limit also applies for the case of a nonlinear source term defined by a power law, as it coincides with the standing solitary wave of the Korteweg-de Vries equation (Regenauer-Lieb et al., 2013b; Veveakis and Regenauer-Lieb, 2015).

III.1.3 Cross-diffusion waveforms

Cross-diffusion waves bear similar characteristic to quasi-solitons encountered in optical systems where chromatic dispersion is strong leading to anomalous dispersion patterns that, unlike solitons, come in discrete portions (Paschotta, 2008).

A classification of the solution of equation (2) for the discussed case of cross-diffusion has been presented recently (Tsyganov and Biktashev, 2014) and a broader definition of quasi-soliton wave has been adopted. Quasi-solitons feature complex dispersion relationships, where the wave velocity of individual waves (phase velocity) have different velocity to their smooth envelope wave groups (group velocity). This leads to complex phenomena similar to envelope solitons in the nonlinear Schrödinger equation shown (Tsyganov and Biktashev, 2014) to be akin to the linear cross-diffusion equation in equation (2). They consequently differ from the classical soliton solution in that the shape and speed of such waves in the established regime does not depend on initial and boundary conditions and is fully determined by the material parameters of the medium that they travel through.

Wave packet solutions emerge when the wave energy is concentrated around a finite wavenumber. When the dispersion is weak, the envelope travels at approximately uniform group velocity. The addition of linear cross-diffusion can create strong non-linear dispersion, where - depending on the coefficients - three different wave types have been classified: (1) fixed-shape propagating waves, (2) envelope waves, (3a) multi-envelope waves, and (3b) intermediate regimes appearing as multi-envelope waves propagating as fixed shape most of the time but undergoing restructuring from time to time (Tsyganov and Biktashev, 2014).

The collisional behaviour of these classes of cross-diffusional waves is complicated. When the cross-diffusion wave hits an interface, or another cross-diffusion wave the amplitude of the quasi-soliton changes and there is a temporary diminution of the amplitude or in extreme cases an annihilation. In most cases they recover gradually their original form. This is another feature that is different to true solitons which do not change on impact. In two-dimensional systems additional complexities arise as they can penetrate or break on collision, reflect into different directions and lead to complex patterns (Tsyganov and Biktashev, 2014) imprinted in the damage zone of seismogenic faults.

Of specific interest for the earthquake problem is the aspect of the fate of the accelerations carried by the waves (Regenauer-Lieb et al., 2019) for cases where the wave reduces in amplitude on collision or annihilates after collision. In experiments carried out with collapsing puffed rice particles (Guillard et al., 2015) cross-diffusion waves were detected (Hu et al., 2019) which release an accoustic emission when annihilating after collision with the top surface. We propose that for a homogeneous plastic zone as shown in Figure (2), the first wave-collisions on the symmetry axis of the future fault plane converts the energy loss of the cross-diffusion wave into local material damage that will act as a seed to future wave interactions and systematically grows the fault. For the case of heterogeneous materials the collision of cross-diffusion waves with internal or bounding material surfaces would cause an alternative seed for the nucleation of earthquakes.

IV Discussion and conclusions

We have introduced a new approach to THMC instabilities using the simple concept of cross-diffusion. The effect of cross-diffusion can, for critical conditions, lead to the formation of cross-diffusional waves which, using the Helmholtz decomposition of the equation of motion, can be decomposed into P- and S- cross-diffusion waves (Regenauer-Lieb et al., 2019). The analysis suggests that cross-diffusion provides a crucial link to allow cross-scale coupling of the multiphysics processes potentially leading to earthquake instabilities. For critical cross-diffusion coefficients and associated reaction terms, an energy cascade of acceleration cross-diffusion pressure waves (diffusion P-waves), from small scale thermo-chemical (TC) and chemo-mechanical (CM) dissipative waves to meso-scale hydro-mechanical (HM) dissipative waves, can be triggered by a geodynamic driving force. This in turn may nucleate shear-cross-diffusion waves in the form of a thermo-mechanical (TM) shear-wave instability.

The multiscale THMC-waves initially are low-energy release volumetric diffusion P-waves free of kinetic energy, but as they trigger thermo-mechanical cross-diffusion shear waves (diffusion S-waves), they can tap into a significant portion of the stored elastic energy around the fault. This process may, for a critical set of reaction-diffusion parameters, ultimately lead to a substantial energy release sufficient to communicate via cross-diffusion with their elastic counterparts (Cartwright et al., 1997), and an earthquake instability occurs.

Numerical solutions of this stiff problem are difficult, and it may be useful to consider weak-formulations by a reduction procedure through adiabatic elimination of the fast cross-diffusion process into an effective cross-diffusion formulation by time integration to the slow self-diffusion time scale (Biktashev and Tsyganov, 2016). Specifically, the upscaled solutions can be regarded as long-wavelength (infinite/quasistatic time-scale limits) of the travelling cross-diffusion wave solutions. The dominant wave number is described by a critical ratio of effective self-diffusion coefficients. The relationship between the long wavelength soliton cnoidal solution of the Korteweg-de Vries equation (Regenauer-Lieb et al., 2013b; Veveakis and Regenauer-Lieb, 2015) and the quasi-soliton solution with oscillating tails can be derived by the Jacobi elliptic functions method of solution. If the modulus of the Jacoby elliptic function asymptotically approaches to zero the solutions shown in Fig. 3 can be obtained (Wang et al., 2018).

These quasistatic spatially inhomogeneous solutions are also known as Turing patterns (Turing, 1952). The upscaled equations do not necessarily honour the co-dependence between the diffusion and reaction rate constants of the cross-diffusion process (Hu et al., 2019; Vanag and Epstein, 2009). Turing patterns appear because around the bifurcation point in the reaction-self-diffusion system, linearly unstable eigenfunctions exist that are growing exponentially with time. In order to regularize the problem, finely tuned non-linear terms need to be introduced in the reaction-diffusion equation or, more conveniently, cross-diffusion-like terms need to be added. We postulate here that THMC-Turing patterns are indeed the multiscale patterns observed in nature (Sethna et al., 2001). They offer themselves as an ideal tool for inversion of the effective self-diffusion THMC coefficients and their implicit reaction rates. Interpreting geological structures in terms of these process parameters will allow identification of principal processes underpinning the earthquake mechanism.

Theoretical considerations on the nature of cross-diffusion waves may help in the design and analysis of laboratory experiments. Quasi-solitons are similar to real solitons in that they can penetrate through each other and reflect from boundaries. Differences are that the amplitudes of the true solitons do not change after impact while the dynamics of quasi-solitons on impact is often naturally seen as a temporary diminution of the amplitude with subsequent gradual recovery (Tsyganov and Biktashev, 2014). Another important difference is that the amplitude and speed of a true soliton depend on initial conditions, while for the quasi-soliton they depend on the material parameters. This property offers new avenues for earthquake physics. We are currently investigating these wave phenomena in controlled laboratory experiments and attempt to use the unique relationship of amplitude and wave speed dependence on material properties as a diagnostic tool for data assimilation.

V Acknowledgments

This work was supported by the Australian Research Council (ARC DP170104550, DP170104557) and the strategic SPF01 fund of UNSW, Sydney.

VI References