Cross-Diffusion Waves as a trigger for multiscale, multiphysics Instabilities: Application to earthquakes

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 the coupling of Thermo-Hydro-Mechano-Chemical reaction-diffusion equations in a THMC diffusion matrix. We show that the off-diagonal cross-diffusivities can give rise to a new class of waves known as cross-diffusion waves. Their unique property is that for critical conditions cross-diffusion waves can funnel wave energy into 5 a quasi-stationary wave focus from large to small-scale. The equivalent extreme event in ocean waves and optical fibres leads to the appearance of ’rogue waves’ and high energy pulses of light in lasers. In the context of hydromechanical coupling, a rogue wave would appear as a sudden fluid pressure spike on the future fault plane. This is here interpreted as a trigger for the ultimate (shear) seismic moment release.

Magnitude Number of Earthquakes Figure 1. Earthquake frequency-magnitude histogram prepared from the Global Instrumental Earthquake Catalogue, Version 7.0 -released on 9/4/2020 by the ISC-GEM (Di Giacomo et al., 2018;Storchak et al., 2015Storchak et al., , 2013. Global quakes have a fractal power law relationship of log(Number of Earthquakes) = 10.23 + 1.06 (Magnitude) with an R 2 = 0.98. The log-log relationship between numbers of earthquakes and their moment magnitude on all sizes of completeness of the catalogue suggests a simple underlying reason why the physics of the very small influences the physics of the very large in a multifractal cascade of instabilities. Sethna et al. (2001) postulates that this can be explained by the existence of a critical thermodynamics force above which an earthquake starts to slip. Accordingly, plate motions self-regulate (Sornette and Pisarenko, 2003) to be exactly at this critical point for potential failure at all scales. systems (Ball, 2012) and their role in thermodynamic far-from-equilibrium systems was originally described by Prigogine 20 and co-workers (Kondepudi and Prigogine, 1998). An application of self-diffusing reaction-diffusion equations to mineralising systems has been proposed recently (Oberst et al., 2018). This contribution investigates geological applications of the possible relation to the reaction-diffusion wave phenomenon due to the new addition of the cross-diffusion term proposed in part 1 (Regenauer-Lieb et al., 2020). Along the same vein, the paper also attempts to give a somewhat simpler description of the theory from a chemical perspective putting part 1 into context with other more familiar theories as well as recent further 25 developments thereof.
The characteristic time and length scales of the thermodynamic THMC processes can be evaluated by calculating the re-100 laxation time of a random perturbation from the THMC reaction term R i . Characteristic time-and length-scales emerge from quasi-steady-state (time-independent) solutions of the diffusive (relaxation) response of the reaction-diffusion system characterised by the equation where C i stands for the diffusing species/processes (e.g. temperature, fluid or mechanical pressure, chemical species) and ζ 105 for the respective diffusivity with the subscript i denoting the THMC processes. Linearising the diffusion equation to first order we obtain 5 https://doi.org/10.5194/se-2020-149 Preprint. Discussion started: 29 September 2020 c Author(s) 2020. CC BY 4.0 License.
The time scale or the relaxation process is now entirely characterised by the self-diffusion constants ζ T,H,M,C and the reaction rate constants R i .

110
2.1 Self-Diffusion length/time scale For a discussion of the length scale encountered through the diffusion of a local reaction source term of the THMC process (i.e. local heat release, fluid or mechanical pressure source term, chemical reaction) we treat the concentration C 0 released by the reaction at t 0 as an initial condition and set R i = 0 for the remainder of the process t > 0 . We use this instantaneous source as a point source in an infinite plane. This linear partial differential equation requires one initial condition (local release 115 of concentration C 0 as a point source due to the reaction at x = 0) and two boundary conditions for solution. The boundary conditions are C i = 0 at plus and minus infinity. Mass (momentum, energy and entropy) balance requires that at any time the concentration (mechanical/fluid pressure, heat content) is: If we assume that the initial and boundary conditions are unchanged by a scale change, as required by the coarse graining 120 assumption, the well-known solution method is obtained by scaling the equation with a characteristic diffusional length scale.
This diffusion length scale defines the characteristic length scale of propagation of the information from a local perturbation, e.g. for chemistry an instantaneous point source in the concentration of the diffusing species being released at the origin x = 0 at t = 0 in a given time t d . For fluid flow and mechanics this would be a local source of fluid or mechanical pressure, and for temperature this perturbation would be a local heat source. Any scale that is significantly larger than the diffusional length 125 scale is considered to be unaffected by the diffusion front for a given diffusion time t d . We will later on reinterpret the diffusion length as an uncertainty measure. The characteristic diffusion length scale is: which is used as part of the scaling method for solving Eq. 2 described in standard textbooks e.g. (Crank, 1975). Using the previously described initial and boundary conditions, the diffusion front spreads radially in a concentration profile The solution shown in Fig. 3 illustrates the fast decrease of the amplitude with increasing time coupled with an increase in L d . The normalisation over 1 L d √ π in Eq. 5 ensures that the area under the curve always remains the same and satisfies Eq. 3 thus defining a Gaussian scale space. The Gaussian probability density function (here called a Gaussian wavelet) is often used in probability theory where the diffusion length is related to the standard deviation, its square to the variance and the function 135 centroid to the mean.

Reaction-Diffusion length/time scales
So far, we have discussed a point source reaction as an initial condition at time t 0 and x = 0 with no additional time-scale other than the self-diffusion process afterwards. This leads to a decaying and broadening diffusion wavelet (Fig. 3). When considering an active source term, the solution turns into a propagating wavefront which was first discovered in 1906 by 140 Robert Luther. An English translation of the original article appeared in the Journal of Chemical Education (Luther, 1987).
The same phenomenon was rediscovered 30 years later and is now known as Fisher or Fisher-Kolmogorov equation (Adomian, 1995). Fisher originally discussed the reaction-diffusion equation to calculate the propagation of a mutant virus in an infinite domain (Fisher, 1937). Following Showalter and Tyson (1987) we recast Fisher's arguments for calculating the minimum speed of the propagation of the mutant gene at long time scales into a discussion of a propagation of a chemical reaction front. For 145 the example discussed in Fig. 3 consider the following autocatalytic chemical reaction: where A is the reactant and C is the catalyst, n is the autocatalysis order with the number of reaction steps with n = 1 being the elementary reaction (Valero and Moyano, 2017). The square bracket indicates the number of moles. The reaction rate is dt with k being the first-order reaction-rate here defined with respect to the reaction product rather than the 150 reactant, hence the positive sign. This reaction depends on the concentration of only one reactant. The system may have other reactants, but these are not influencing the rate of the reaction. In chemistry this situation is known as a first-order reaction.
As we are interested in a moving reaction front, we may consider just the initial response of the reaction. Thus instead of the typical sigmoid solution, which considers the depletion of [A] in the course of the reaction, we may consider a non-physical infinite pool, where [A] is not depleted by the reaction, and integrate the first-order reaction by an exponential growth of the 155 reaction product, with The minimum velocity of the propagating reaction-diffusion wavelet that is triggered by an initial delta-function reaction at time t 0 can be evaluated by extending the solution in Eq. 5 to consider the exponential reaction product and obtain: 160 Separating variables , and substituting the separation into Eq. 8 we obtain an ordinary differential equation (ode) 1 @0s @2s @4s @6s @8s @10s x Figure 4. Illustration of the propagating wavefront of the Fisher-Kolmogorov equation which is Eq. 6 with a nonlinear source term Ri = kiCi(1 − Ci). Kolmogorov et al. (1937) showed that any initial concentration vanishes for large x and evolves to a travelling wave solution with a minimal velocity v = √ 4ζiki . The reaction-diffusion exhibits a bistable equilibrium. One stable equilibrium is achieved when the initial concentration is below the activation of a wave and the system is resting and the other is the travelling-wave solution with a minimal velocity v. We show here the non-dimensional solution for the example of a tanh initial condition (black curve) which quickly converges to a travelling waveform with a non-dimensional speed of 2. Various solution techniques exist. A convenient method is to turn the partial differential equation (Eq. 2) into an ordinary differential equation by using Chebyshev polynomials (Towers and Jovanoski, 2008). where ζ is the diffusivity. This shows that for large t and large x the last two terms of Eq. 10 vanish, and a constant wavefront which is a lower bound of the steady front velocity at large t. We can use this speed as a quasi-steady-state solution of the composition wave triggered by the reaction and obtain a quasi-steady state diffusion length scale of L d = ζ v that propagates with the wavefront. This diffusion front is the only information, that is carried by the wave at long length and time scale, and 170 highlights one of the fundamental differences between reaction-diffusion waves and elastic waves. Reaction-diffusion waves quickly diffuse information of the initial condition and at long time scales are only characterised by the competition between reaction rates and diffusion rates. Pure elastic waves, when propagating without damping, diffraction and scattering, carry the full source information through their entire travel path.
The above example of a chemical autocatalytic reaction illustrates the key features of reaction-diffusion equations in many 175 scientific disciplines and applies to all C i 's and R i 's of THMC systems. Adding an active non-linear source term (e.g. Fisher-Kolmogorov equation) into the diffusion equation can lead to the interesting phenomenon of the generation of a self-oscillatory excitation wave, where after Fisher's work on the topic progress was mainly made in the Russian literature triggered by the seminal work of Kolmogorov et al. (1937) discussed in Fig.4.
In Russian literature, the term "autowave" was introduced (Ostrovskiȋ, 2015). The autowave phenomenon is well-studied for 180 understanding the electrical nerve impulse that leads to contraction of the heart muscle (Antonioletti et al., 2017). However, it constitutes a fundamental class of waves encountered in all reaction-diffusion systems in physics, biology, and chemistry (Vasil'ev, 1979). The principal difference to classical wave equations that are based on hyperbolic differential equations, is that the autowave phenomenon arises from a non-linear source term in parabolic equations. When adding a cross-diffusion term autowaves can turn into standing waves (Berenstein and Beta, 2012) having completely different properties to classical standing 185 wave solitons. This difference defines a new class of dissipative waves which, according to Tsyganov and Biktashev (2014), presents an entirely different "world" to the waves encountered in integrable conservative systems such as the elastic-wave phenomenon. This phenomenon will be discussed further in the section on cross-diffusion waveforms where we will discuss the most important property for nucleation of earthquakes. This is the potential of multiscale funnelling of wave energy from the environment into a localised standing wave. Before going there, we continue with a discussion of the reaction-diffusion 190 equation without the cross-diffusion term in order to systematically describe its basic characteristics.
The reaction-diffusion time-scales of a single reaction-diffusion equation can be expressed by a generalization of the Damköhler number for chemical reaction-diffusion processes to all THMC couplings. The Damköhler number describes the ratio of the diffusion over the reaction time or the equivalent ratio of the reaction over the considered diffusion rate.

195
This ratio now defines a new time-scale and replaces the purely diffusive time-scale discussed in the spreading wavelet shown in Fig. 3 with a propagating wavefront solution that self-supports its shape. This is illustrated here only for an infinite autocatalytic source term in the reaction-diffusion system (Fig. 4). For a finite autocatlytic reaction the reactive source term has a growth function which is sigmoid, also called a logistic function. This means that the reactant initially grows exponentially, similar to the infinite source, followed by linear reactant growth and a final zero growth branch. These general solutions imply different 200 shapes of the self-supporting propagating wave and, depending on the value of the diffusivities, also a finite life-time of the wave.
This geologically more relevant situation is described in Molotkov and Vakulenko (1993). The authors describe generalised reaction-diffusion systems and find that the wave behaviour depends on only three parameters. For small wavefront curvature the autowave is described by the two parameters of the infinite source term solution discussed above, while for the more general 205 case the normal velocity of the front may contain the front curvature as an additional parameter. The additional dependence of the wave function on the curvature of the concentration field arises because of the fact that ∂Ci ∂t is proportional to the curvature of the wavefront. For regions of the wavefront where the curvature is negative the concentration must decrease at a rate proportional to the magnitude of the curvature. Conversely, the concentration must increase where the wavefront curvature is positive.

210
Summarizing the above findings, we can now characterise the reaction-diffusion thermodynamic system by five key features: i.) a bistable or multistable (for several reactions) region with a stable stationary mode and a mode for the nucleation of propagating autowaves above a critical activation threshold; ii) in the activated state, the wavefront separates two regions, a local region characterised by the particular THMC diffusional length scale L d affected by the reactions R i , and a large region at > L d which is outside the reaction-diffusion wave; iii) for long time scales, the wavefield is governed by characteristic 215 self-oscillatory motions which for bistable systems are described by just three parameters. For multistable systems chaotic oscillations are expected (Molotkov and Vakulenko, 1993). For the analysis of this complicated system we will propose to use perturbation theory and illustrate the approach through some basic concepts of signal processing. In the bistable system the Fisher-Kolmogorov wave is a self-propagation dissipative wave at a characteristic wave speed whose lower limit can be quantified by the square root of the Damköhler number times the diffusivity normalised by the characteristic diffusion length 220 scale L d ; (iv) the propagating wave exponentially scatters information about its initial condition, and the wavefield only carries information about the dissipative properties into the far-field; this is an important differentiation to waves in the conservative system (e.g. elastic waves) where the wave at a long-distance still carries information about its initial conditions; (v) the wave speed thus becomes a fundamental material constant defined by the rates of the dissipative THMC processes as: 225 As the spreading wavefronts are self-supporting and can propagate upwards in scale, we propose that this material velocity not only applies to the above discussed chemical reaction-diffusion equations but to all reaction-diffusion equations of the THMCcoupled system. In this proposition, the propagating multiscale and multiphysics waves provide the capacity to link the different THMC-length-scales and could explain the multifractal nature of earthquakes. The approach allows a significant simplification of the earthquake physics problem as the exponential rate of approaching the Kolmogorov limit of a self-oscillating wavefront 230 shown in Fig. 4 can be simplified by replacing the nonlinear source term through a set of at least two linear coupled partial differential THMC diffusion equations which will be discussed in the next section. At a given scale, the Kolmogorov wave velocity limit allows a characterisation of the important physics in terms of the wave velocity. With the autowave approach, one can turn any non-linear perturbation of the local source into a linear propagating waveform only governed by the dissipative material properties. Autowaves will by themselves recover a characteristic wavefield dictated by the reaction-diffusion rate 235 constants.
This characteristic behaviour is used, for instance, in medicine where autowaves are encountered in many fields. The electrical nerve impulses that drive a regular heart beat are an example (Antonioletti et al., 2017). The authors describe how the characteristic recovery of the autowave waveform after a random perturbation can be used, for instance, for defibrillation strategies (a small electrical stimulus applied through a pacemaker) for treatment of life-threatening heart arrhytmia. In this sense, 240 earthquake physics might profit from an understanding of the partial differential equations developed in mathematical biology, epidimiology (wave-like propagation of viruses) and other biomedical applications for which numerical open-source tools are available. One such tool is the "heart beat box" (Antonioletti et al., 2017) which uses a mathematical formulation of the human heart in terms of a coupled electro-mechanical reaction-diffusion equations similar to the coupled reaction-diffusion equations discussed above.

245
The foregoing results were based on the generalisation of a first-order autocatalytic reaction. Very few practical examples will be first order, and we need to consider the generalised case of a reaction that depends on the concentrations of a second-order reactant or any number of additional reactions and their equivalent THMC processes. To extend the discussion to more than one reaction, we apply the Fisher reaction term to more than one interdependent reactant. That means in chemical terms that the concentration of second-or higher-order reactions depends on the concentration of a second-or higher number of species. This 250 will add additional dynamics to the system response. To isolate the effect of this additional feedback it is useful to neglect the spatial response and neglect diffusion. Mathematically, such an interdependent coupling of reaction rates is expressed by just considering the reaction term for two cross-linked autocatalytic equations. This is known as the Lotka-Volterra predator-prey model (Lotka, 1920), which beautifully illustrates the generic behaviour of a coupled system in time. The interaction with the spatial response will be discussed later.

Periodicity in time: Lotka-Volterra Waves
The Lotka-Volterra predator-prey model couples two autocatalytic reactions A + X → 2X prey (X) have an infinite supply of food (A) and multiply Y + X → 2Y predators (Y) eat prey (X) and multiply as a consequence .
Coupling between the two autocatalytic reactions occurs because the predators can only multiply if they eat prey. Similar to the above assumption for the derivation of the limiting wave speed, we assume that the food supply A for the prey is infinite 260 and without predators there would be an exponential growth of prey. In chemical terms for the molar concentrations, this leads to the following rate equations (Lotka, 1920).
where k 1 and k 2 are reaction rates of the autocatalytic reactions derived from empirical/phenomenological postulates. To close the system a third rate k 3 must be introduced that quantifies the rate of death of predators as the death rate of prey is 265 already included in the reproduction equation for the predator. We can evaluate the equilibrium points of the system by setting dt = 0 and obtain following condition for the rates: Another equilibrium point is where the predators consume all prey leading to the extinction of both species. Integrating the Lotka-Volterra ode numerically reveals the fundamental oscillatory behaviour of coupled reactive systems shown in Fig. 5.

270
The Lotka-Volterra oscillator applies to second-order and higher-order reactions where multiple oscillators may be encountered. We posit here that the fundamental analysis can be transferred from the classical biological system (predator-prey, infectious diseases, etc.) and chemical reactions to all coupled reactive source terms of the THMC reaction-diffusion system.
When looking for a possible Lotka-Volterra oscillator in earthquake dynamics, one would start with trying to identify first a perfectly regular oscillator which would involve only two coupled THMC equations. To then generalise the approach one 275 would continue to the generalised Lotka-Volterra oscillator, such as the three-species oscillators, which explain the transition to chaos. An example is shown in Fig. 12.4 in Flake (1998).

Lotka-Volterra type earthquake sequences
There are examples of characteristic earthquakes with some regularity in their recurrence, e.g. the Parkfield example (Wiemer and Wyss, 1997)), and also the Episodic Tremor and Slip (ETS) events recorded in Japan, Cascadia and Hikurangi subduction 280 system (Gomberg, 2010). These have been interpreted with a perfectly periodic oscillator model (Poulet et al., 2014b) shown in Fig. 6. It was shown that a coupling of the temperature reaction-diffusion equation with the pressure reaction-diffusion equation Veveakis et al., 2014) can lead to the characteristic perfectly periodic Lotka-Volterra-type oscillatory response in the temperature-pressure plot (Fig. 6). The authors also show an example where, by considering an additional oscillator, a transition to chaos can be modelled (Poulet et al., 2014b). The basic element of the model is a chemical decomposition 285 function of the type AB A + B. This dissolution decomposition reaction is different to the autocatalytic reaction described above. The equivalent autocatlytic element in this decomposition reaction is introduced through shear heating feedback during Figure 6. The characteristic period of the Cascadia ETS sequence interpreted with a reaction-diffusion oscillator model (Poulet et al., 2014b).
(a) The phase diagram between A and B corresponds to the slow creep of the serpentinite slowly raising the temperature to the critical level for the onset of the dehydration reaction. The segment between B and C corresponds to the dehydration reaction that coincides with the tremor and slips event. The following segment between C and A is the diffusion-dominated field. (b) shows the corresponding normalised strain rate with the sharp peak corresponding to the tremor event.
fault slip, triggering thermally induced dissolution. Poulet et al. (2014b) numerically investigated the potential candidate source mechanism of a thermally induced dissolution reaction of serpentinite in the fault zone originally proposed by Obara (2002).
In that model solid serpentinite (phase AB) decomposes into a solid phase A (Antigorite) and a fluid phase B (water) upon 290 application of a heat source from shear heating.

Periodicity in space: Cross-diffusion Waves
It is encouraging that a simple chemical reaction model of the normally fine-grained serpentine crystal at say millimetre scale can be used to model an ETS instability at plate scale, say around 100-1000 km scale. Communication of information over eight to nine orders of magnitude may hence be possible under special circumstances. This leads to the proposition that progress 295 can be made by investigating the multiscale physics of the THMC reaction-diffusion system in the Earth in further detail. In addition to the timing and displacement information obtained from GPS and seismic stations, the spatial information related to the diffusion term could be investigated to verify the model. A key observable would be the propagating wavefront that would be expected from the above discussion on reaction-diffusion time scales.
In the case of the serpentinite decomposition reaction, the diffusive wavefront would propagate normal to the main fault plane 300 as the chemical decomposition reaction involves significant volume reduction which would lead to a contraction of the central fault plane and a characteristic width of the fault plane. This process zone is not exposed in active subduction systems but can be found in geological exposures of fossil faults. An excellent example is the Glarus Thrust which has been modelled by the Lotka-Volterra-type oscillator approach using a carbonate decomposition reaction (Poulet et al., 2014a). The model includes diffusion and therefore the expected spatial response can be tested as well. The model is shown to be capable to reproduce the 305 process zone around the central fault plane which is documented in many fault zones (Chester et al., 2013). Another textbook example is the Punchball Fault in California where a series of deformation mechanisms were described (Schulz and Evans, 2000).
Synchronous or staggered deformation mechanisms are a general observation (Vermilye and Scholz, 1998) and require expansion of the chemical decomposition reaction model to include other important meso-to macroscopic deformation mech-310 anisms. Additionally, microstructure, geometric complexity, and multiple deformation mechanisms lead to a stochastic element at the meso-scale whose role in the transfer of the information from micro-to large-scale in real geological applications is not yet discussed. In the following we will discuss whether these meso-scale stochastic processes lead to a loss or an enhanced coupling from small to large scale.
The problem can be addressed by a statistical mechanics approach to the coarse-graining problem. In Fig. 2 we proposed 315 that the missing element for considering mesoscopic complexity in a physics-based earthquake model is a meso-scale approach that captures the link between the vast differences in the diffusional length scale of THMC processes shown in Table 1. In part 1 (Regenauer-Lieb et al., 2020) we proposed that a meso-scale approach can be developed by decomposing the large-scale reaction term R i into a meso-scale reaction term which requires a meso-scale cross-diffusion term for mass (momentum, energy) balance. This is because the aforementioned cross-scale coupling introduces meso-scale source/sink terms in the individual 320 conservation laws identified by r T , r H , r M and r C , respectively. The conservation laws must, therefore, be extended to close the equations and allow a multitude of THMC processes to occur simultaneously, which introduces cross-diffusion fluxes.
Here we derive the cross-diffusion approach from a different angle by following the coarse-graining operation. Techniques for the quantification of the uncertainty reduction through the coarse-graining operation are the main subject of statistical mechanics (Sethna, 2006). The basic principles can be discussed based on Shannon's entropy, Heisenberg uncertainty relationship, 325 renormalisation theory, and the fluctuation theorem (Evans and Searle, 2002). For a more in-depth discussion and a broader perspective on the different techniques we recommend textbooks providing distillations of 50 years of statistical mechanics (Sethna, 2006) and material science applications (Balluffi et al., 2005). Here we present a functional analysis approach of probability theory applied to order parameter fluctuations. The first part of the analysis is often used in signal processing and image analysis (Buades et al., 2005). The following discussion, therefore, explains the convolution filter interpretation of the 330 cross-diffusion waves discussed in part 1 (Regenauer-Lieb et al., 2020).

Gaussian Wavelet convolution as a blurring filter
In terms of uncertainty quantification the diffusion wavelet in Eq. 2 and shown in Fig. 7 is identified as a univariate Gaussian probability density function, also called a Gaussian kernel in image processing (Buades et al., 2005). Due to the linear nature of the partial differential equation we note that the following analysis is valid for any number of interactions. Fig. 7 shows an 335 example of interaction of two diffusion wavelets with the probability distribution functions having an opposite sign of velocities Figure 7. Consider a thermodynamic system as a planar surface from a macroscopic view. At the mesoscopic view, it is made up of subvolumes or more accurately subplanes in the current mathematical formulation. A mesoscopic subvolume is highlighted where two propagating diffusion wavelets interact. The discrete system represents the microscopic view where red and blue dots represent individual discrete chemical molecules or predators and prey or indeed generalised thermodynamic 'micro-engines' in an activated oscillatory quasi-steady state. We are interested in the assessment of the effect of the mesoscopic collision of micro-engines on the macroscopic system. From an equilibrium statistical mechanics perspective, we are therefore not interested in the local dynamics of molecules/micro-engines but use a convolution operation to describe the mesoscopic effect. This approach allows investigation of the meso-scale dynamics which feeds into the macroscale continuum mechanics approach where sufficient time has elapsed so that the flux Q becomes constant and the system has reached its maximum entropy production state (Collins and Houlsby, 1997). The system can then be described by the theory of thermomechanics, ie.
the theory of thermodynamics with internal variables (Maugin and Muschik, 1999;Jacquey and Regenauer-Lieb, 2020). This approach also yields the familiar dissipative material properties/transport coefficients of the classical continuum approach. For constant thermodynamic force flux pairs we obtain for instance: the thermal conductivity of the material, for constant temperature/heat flux (T), permeability for a constant pressure/fluid flow pair (H), viscosity for a constant velocity/momentum flux pair (M), and the chemical diffusivity for a constant chemical potential/diffusion flow (C) pair (Oettinger, 2005).
implying a collision of the two wavelets due to their opposing direction of travel (Eq. 13).
where µ f and µ g are the function centroids of the the function f (x) and g(x) known as the mean in probability theory. Since both wavelets are moving towards each other we would like to use a Lagrangian reference frame, here arbitrarily chosen to be 340 f (x − τ ), to assess how the shape of this function is modified by the passage of the function g(τ ). The convolution operation is commutative, and we can also use the opposite reference frame. The convolution operation is defined by: 16 https://doi.org/10.5194/se-2020-149 Preprint. Discussion started: 29 September 2020 c Author(s) 2020. CC BY 4.0 License.
where τ denotes a translation in the positive x-direction. Due to the translational symmetry of the diffusion wavelet, it is convenient to use the Fourier transform F as the eigenfunctions of the translation τ operation are complex plane waves.
where k is the wave-vector of the Fourier transform. The wave-number and wave-vector are defined as k = |k| = 2π/λ. The convolution operation of function f (x) and g(x) in Fourier space is the product of the Fourier transforms F (f (x))F (g(x)).
Applying the same Fourier transform as in Eq. 19 to F (g(x)) and performing the product with F (f (x)) we obtain: 350 The first exponential term is the mean amplitude of each Fourier wave-vector which can be decomposed as a real sinus-wave and imaginary cosinus-wave that contains the phase information. The second term is a Gaussian probability density function of the wave-vector.
This result is very useful as it can be used to analyse the behaviour upon the collision of two diffusion wavelets. The first conclusion is that the convolution results in another Gaussian wavelet as can be seen from the comparison of Eqns. 20 and 19.

355
The second conclusion is that the convolved amplitude is reduced by the sum of the function centroids µ f + µ g in the positive x-direction and the breadth of the wavelet is enlarged by the variance (L d(f ) + L d(g) ) 2 of the probability distribution function.
The result that the convolution operation reproduces a similar function to the original function is called self-similarity. The Gaussian function is self-similar and defines a linear Gaussian scale-space where space-time relationships can be normalised by the diffusion length. This is an important aspect of coupling across the scale. Another aspect of the Gaussian wavelet is 360 that the 1-D approach presented here can be extended and superposed to any dimension as the function is separable, i.e. it can be applied to any dimension independently. In principle the technique can be extended to anisotropic diffusion, where the introduction of a non-linear space varying Gaussian filter may lead to additional problems discussed next (Kichenassamy, 1997).

Gaussian Wavelet convolution as a sharpening filter
365 So far we have discussed the situation where the meso-scale convolution of the two diffusion wavelets appears as a smoothing filter of colliding Gaussian wavelets, where high frequencies are quickly reduced and the high-frequency wavelet gets preferentially diffused. This class of Gaussian low pass filter problems is not of interest for the nucleation of earthquakes as the small scale local instabilities triggered by chemical reactions are filtered out. The opposite situation of a sharpening filter is, however, very interesting. A sharpening filter is an expression used in image processing when applying a deconvolution operation using 370 a Dirac delta convolution kernel minus a Gaussian blur kernel. This is effectively an anti-diffusion filter, also known as an uphill diffusion operation, which may lead to problems as it may violate thermodynamics without applying further restrictions (Ulmer, 2010). This problem is well known in image processing for modelling anisotropic diffusion where it is known as the Perona-Malik paradox (Kichenassamy, 1997).
However, the sharpening filter problem can be made thermodynamically consistent. The method has been described in part radiate energy into the environment and regularise the local inconsistencies through their interaction with a stabilising constant large-scale thermodynamic flux. The communication mechanism between the meso-scale and macro-scale is a travelling acceleration wave. We have shown through mass and momentum balance that a local reaction term requests a local cross-diffusion term. The approach in part 1 (Regenauer-Lieb et al., 2020) was based on a physics and mechanics of solids approach which es-380 sentially looks at the meso-scale negative diffusion process from a macroscale continuum mechanics perspective characterised by a positive diffusivity. Can we come to the same conclusion by looking at the meso-scale from a chemical-scale perspective by considering the mobility of atoms underpinning the meso-scale process in the light of the global changes of the free energy?
The observation of spontaneous unmixing from a fully diffused, intermingled, and unstable thermodynamic phase of two different phases, e.g. solid and fluid (spinodal decomposition), is mathematically equivalent to an uphill diffusion problem.

385
When interpreting the phenomenon through a linearised Fickian diffusion Ansatz (Eq. 2) we come to the perplexing conclusion that Eq. 2 has a thermodynamically inconsistent negative diffusivity (uphill diffusion). Note that the driving force for spinodal decomposition is not a chemical reaction but a reduction in the bulk free energy of the homogeneous mixture, which is thermodynamically unstable. The breakthrough in understanding this process in terms of a thermodynamically consistent approach to the apparent uphill diffusion phenomenon of spinodal decomposition is attributed to Cahn (1961) and coworkers.

390
A detailed description of the approach is given in Balluffi et al. (2005). The solution is to consider a meso-scale length scale at the spinodal reaction front (hereafter called the spinodal) and express the diffusion process as a product of the diffusive mobility of atoms M , which must be a positive term to satisfy the second law and a gradient term of the large-scale generalised chemical potential. The gradient in the generalised chemical potential, in turn, is made up of a product of a meso-scale negative second-order derivative of the free energy f of the mixture concerning the concentration C and a positive gradient of the 395 concentration. This situation leads to a non-linear interdiffusivityD = M ∂ 2 f ∂C 2 which is negative at the spinodal.
There are two ways to approach this thermodynamic problem which can be solved by considering the meso-scale thermodynamic couplings. The first approach is similar to the one presented in part 1 (Regenauer-Lieb et al., 2020) but developed from a chemical perspective. For this, we decompose the nonlinear interdiffusion coefficientD into n-coupled interdiffusion 400 coefficientsD ij expressed as a product of two matrices L ik and Υ kj (Balluffi et al., 2005), where L ik are the classical macro-scale chemical diffusivities, also known as the 'Onsager coefficients', and Υ kj are the meso-scale thermodynamic factors that couple local chemical potentials to concentrations. For the specific case where the local, meso-scale coupling coefficients are of mechanical origin, due to the effect of the volumetric strain triggered by a concentration 405 change for instance, Υ kj are known as the 'Vegard coefficients' (Balluffi et al., 2005).
For the least complicated multicomponent system the coupled interdiffusion matrix is We will use the compact notation for the extended general multicomponent case and introduce the stability criterion for the interdiffusion matrix later. For the spinodal decomposition case an instability is triggered when the meso-scale couplings 410 have a negative interdiffusivityD ij and a local spinodal wavefront develops. A negative interdiffusivity leads to an exponential growth of any perturbation exactly like in the previously discussed reaction-diffusion case.
For discussing the criterion for instability we follow the solution strategy originally suggested by Cahn (1961). For this we stay in the Fourier space and consider the general linearised diffusion problem:

415
whereby C m is the average composition used to lineariseD into the components of the two phases. A(k) are the amplitudes of each Fourier mode k being Inserting the chosen interdiffusion couplings into the generalised solution Ansatz we can estimate the response of the system through a perturbation analysis of the wave-vector in Fourier space whereby and R(k) is the amplification factor. This generic perturbation analysis has the following solution: This analysis divides the spinodal decomposition model into two areas. For the local scale around the spinodal composition wave we expect positive values of R(k) and the convolution of the two inter-diffusing phases is unstable and results in ex-425 ponential growth of composition waves defined in Eq. 24, i.e. a sharpening filter with the maximum wave-vector ∂R(k) ∂k = 0. Away from the spinodal composition waves, the amplification factor is smaller than zero (R(k) < 0), and the system is stable. Figure 8 shows three-time steps in the temporal evolution of a simple Cahn-Hilliard system illustrating dominant wavelength growth at the maximum amplification factor R(k).
To visualize the spinodal wavelength changes in the system, we plot the normalized species concentration in the x−direction.

430
The dynamics increases the spinodal wavelength as the system goes to the steady-state which drives the phase-separation process. For Fig. 8 we solve the Cahn-Hilliard equation using PetIGA in a 2D configuration with 64 finite elements in each direction. PetIGA is an open-source framework for high-performance computing that solves efficiently phase-field equations (Dalcin et al., 2016). The material parameters for the Cahn-Hilliard system allow for phase separation together with spinodal decomposition. Moreover, we use the isometrical analysis method to successfully discretize and solve the equation in its primal 435 form. Such a method allows for high-order, highly-continuous basis functions (Cottrell et al., 2009).
For the generalisation to a multi-component system with N -coupled components, we evaluate the bi-directional response of each coupled system through a complex-valued wave function perturbation. This is also known as an order parameter perturbation analysis cross-diffusion of u: This bi-directional relationship between u-and v-wave, more generally known as the Kramers-Kronig relationship has provided a robust approach to investigate order parameters, phase transitions, and fluctuations (Balluffi et al., 2005;Sethna, 2006).
An important conclusion from this formulation is gained by a comparison of the partial differential Equations 29 and 15. The cross-diffusion problem leads at a quasi-steady state to a periodicity in space in the same way as the Lotka-Volterra oscillator model leads to a periodicity in time at quasi-steady state. The partial differential equations for cross-couplings are the same except that the cross-coupling coefficients in the Lotka-Volterra oscillator (Eq. 15) are reaction rates and the cross-coupling coefficients in the cross-diffusion formulation (Eq. 29) are diffusivities.

Cross-diffusion and its role in coupling of instabilities across scales 455
After discussion of the isolated spatial (cross-diffusion) and temporal (autocatalytic reaction) oscillators we can now extend the above approach for full space-time coupling for multiple reaction-diffusion equations and re-introduce reactions. An exemplary convolution of the mesoscopic oscillators is the self-similar Gaussian wavelet in space and time, which is also the Green's function of the local equilibrium statistical mechanics approach. In the linear case, any number of wavelets can be superposed, and a Fourier transform is an ideal procedure to capture the effect of multiple oscillators and local heterogeneities on the macro-460 system. In part 1 (Regenauer-Lieb et al., 2020) we proposed to generalise the cross-diffusion approach, which is well known in chemistry (Manning, 1970), to more general THMC terms defining cross-diffusion as 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

465
The species j is identified as the cross-diffusion species other than the species i . Introducing a fully populated (N × N ) diffusion matrixD ij , equation (30) can also be written as whereby the classical (self-)diffusive length scale of each THMC process is defined by 4D ii t. The cross-diffusion approach introduces also a new meso-scale coupling length scale that can provide a link between the large scale self-diffusion length 470 scales of the THMC processes. In part 1 (Regenauer-Lieb et al., 2020) we have shown that the cross-diffusion length scale 4D ij t is is defined by the kinetic material properties which can be evaluated by measuring the velocity of cross-diffusion waves in analogy to the generalisation of Eq. 11 for the Fisher-Kolmogorov wave. This leads to the equation for the wavespeed in Eq. 13 defined by the generalised Damköhler number and the cross-diffusion length scale 4D ij t. Using the alternative mechanical formulation the cross-diffusional wave speed turns out to be (Coleman and Gurtin, 1965) 475 where ρ is the density and C is the 4 th -order material stiffness tensor also known as the acoustic tensor discussed in part 1 Any local thermodynamic incompatibility that leads to complex eigenvalues ofD ij must radiate energy in oscillatory instabilities by so-called acceleration waves (Regenauer-Lieb et al., 2020), relaxing to the equilibrium state (Vanag and Epstein,485 2009) to recover the second law of thermodynamics at large scale. 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. The important element of cross-diffusion waves for earthquake physics is their capability to link one thermodynamic force 495 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 before earthquake instabilities. Before discussing the potential earthquake application, we summarise observations from the laboratory and the field. field applications, the processes are often occluded to direct observations, and the self-and cross-diffusion coefficients cannot be easily derived. The verification or falsification of the reaction-diffusion "wave mechanics" approach is facing the following 505 difficulties: (i.) Field observations offer a frozen-in snapshot of the dynamic process and may enable a direct identification of cross-diffusion length scales through observation of the microstructure. However, the assumption is that the reaction-diffusion wave has been given sufficient time to reach its maximum wavenumber as shown in Step 3 in Fig. 8. (ii.) Laboratory measurements offer insight into the dynamics, but they have the drawback that the microstructural processes at the meso-scale are difficult to monitor directly (Schrank et al., 2020), and what is recorded is often an average response. This averaging require-510 ment some heuristic assumptions and simplifications. (iii.) The third and perhaps most important abstraction is that in the laboratory or field reference frame we most often see laboratory or earthquake instabilities as mechanical and not as a chemical, fluid, or thermal instability. The third problem can therefore only be addressed by exploring the solution space mathematically. This paper proposes that earthquakes are a THMC convolution of all of these instabilities that cause volumetric and shear strains due 515 to their different micromechanics. In the following, we will address the three problems sequentially using field and laboratory examples before exploring simple analytical solution applied to earthquake physics.

Diffusion waves at chemical scale
At first sight, we may expect that field examples offer the most direct access to the verification of the cross-diffusion wave hypothesis, however, this is curtailed by following difficulties. The foregoing discussion on the Lotka-Volterra oscillator model 520 compared to the spinodal composition wave highlights the fact that oscillatory responses can either stem from an autocatalytic reaction (Lotka-Volterra, oscillation in time) or an effective uphill diffusivity (oscillation in space) or both. In a geological field example we can investigate the frozen-in spatial response but cannot verify/falsify the theoretical expectations for the time evolution and hence, we cannot quantify the composite effect of the different wave speeds. From a mathematical point of view we may expect two end-member cases. The first case is the temporal oscillator-case that is governing the system as in 525 the Lotka-Volterra model (Eq. 15) and the other is the spatial oscillator-case that is governing the system in the form of waves fixed in space. This spatial oscillator is obtained by swapping the rate constants of the Lotka-Volterra cross-coupling rates by cross-diffusion coefficients as shown in Eq. 29. These two end-member cases are exactly what has been proposed to explain periodic patterns in the field (L'Heureux, 2013).
Accordingly, time-periodic chemical systems have been attributed to a mechanism discovered in reaction-diffusion ex-530 periments with hydrogels where colloidal precipitation patterns with a self-propagating rhythmic growth signature have been recorded (Liesegang, 1906). This special class of Turing precipitation patterns, encountered in reaction-diffusion systems, have been attributed to a large class of oscillatory pattern formation observed in chemistry, physics, biology (shells), medicine (gallstone, cysts, tumours, inflamed tissues) and geology (Nabika et al., 2020). In geology these oscillatory precipitation patterns are encountered in a wide range of geological settings (igneous, sedimentary, hydrothermal and metamorphic) and classified 535 as "Liesegang" patterns (L'Heureux, 2013). Rhythmic Liesegang patterns come in many guises and for widely different rock compositions. Figure 9 shows a typical set of Liesegang patterns encountered in a porous sandstone on the East Coast of Aus-!" #" Although Liesegang patterns have been encountered in many guises and more than 100 years of research on the subject has supplied a wealth of observations, a recent review (Nabika et al., 2020) reveals that the basic mechanism is still controversially discussed. Accordingly, the basic mechanism is attributed to two competing models classified by Nabika et al. In contrast the post-nucleation model follows as a quasi-static pattern which emerges as a long time scale pattern after the reaction-diffusion front has propagated through the medium. The latter model has been used to explain the origin and mechanism for banded iron formations (Wang et al., 2009) and striped zoning in agate (Wang and Merino, 1995). A more recent 550 study of mineralising systems interpreted as long-time scale solution of the Gray-Scott reaction-diffusion system is found in Oberst et al. (2018). The Gray-Scott model neglects the meso-scale cross-diffusion term but considers self-replicating patterns of the mineral system reproduced by inclusion of a non-linear source term. This approach therefore jumps straight to the longtime scale solution as for the equivalent mechanical case of the Korteweg deVries -or the non-linear Schrödinger equations discussed further below.

555
A universal mechanism for Liesegang patterns probably encompasses both end members in a mixture (Nabika et al., 2020).
This conclusion is supported by similarity of the Liesegang patterns to the spinodal decomposition models. In these models the mathematical analysis of the Cahn-Hilliard spinodal composition wave reveals a surprisingly rich field of linear and non-linear wavefronts with possible cross-overs from temporal to spatial wavenumber selection and wave speeds (Scheel, 2017). Such a complex behaviour may be expected from the above discussion that includes the partial differential equation of a  Volterra oscillator in time combined with the cross-diffusion oscillator in space. Due to the complexity of the behaviour we therefore recommend to use the perturbation theory in terms of a Kramers-Kronig formalism (Eq. 28) for a full exploration of the solution space. For direct analysis of sparse data provided by field observations such a complete analysis is, however, not appropriate and simplifications need to be introduced.
In order to do so we need to perform laboratory experiments with analogue materials to gain insight into the timescales. For 565 such experiments, we may switch from a chemical perspective with fluxes of species to a mechanical viewpoint represented by strain and strain-rates. Both can be coupled through an equation of state approach and mixture theory as discussed in part 1 (Regenauer-Lieb et al., 2020). Another approach is to jump from the meso-scale simulation directly to the continuum scale through the introduction of an internal variable and the assumption of local equilibrium. This approach is also known as or thermodynamics of internal variables (Maugin and Muschik, 1999;Jacquey and Regenauer-Lieb, 2020). We will discuss both 570 approaches in the next section.

Diffusion waves at mechanical scale
Oscillatory deformation bands are well documented in deformation experiments of plastic materials such as steel. Perhaps the most well-studied effect is the development of characteristic Portevin-Le Chatelier (PLC) bands which show rhythmic bands of volumetric strain in metals deformed under tensile conditions. The banding is attributed to a damped runaway effect induced 575 by a critical (negative) strain-rate sensitivity (Zaiser and Hähner, 1997). In PLC the transition from smooth to oscillatory deformation is understood as a critical point phenomenon where mesoscopic fluctuations manifest themselves at the macroscale due to strain-rate softening.
A comprehensive review of theoretical approaches to model the phenomenon of an oscillatory material response of plastic materials with a serrated stress-strain response (or jerky flow) can be found in Zaiser and Hähner (1997). According to the 580 review models first have been developed based on a phenomenological approach proposing that the effective stress applied at the boundary of the deformed sample is a non-linear function of the strain rateε which can become negative. In the following, we will interpret the approach at the mechanical scale in the light of the reaction-diffusion approach and THMC coupling concerning the earthquake application.

Oscillatory deformation bands in plastic deformation of materials with internal structure 585
The phenomenological approach where the effective stress weakens with increasing strain-rate is ill-posed as it leads to a runaway effect when a small local perturbation in strain-rate grows without bounds in an ever-increasing manner. Various approaches have been proposed to address this problem. In the context of relating the phenomenological approach to a physicsbased model, the most interesting approach is the one by Estrin and Kubin (1995). The authors proposed to regularise the illposed localised runaway strain by introducing a pseudo-diffusivity; the function that was proposed as an extension of Orowan's 590 equation iṡ where b is the Burgers vector, ρ m the dislocation density, and v d the dislocation velocity and D a diffusivity of the strain ε. This heuristic approach was later shown to be resolved as the 1-D solution obtained from the coarse-graining of the time-evolution gradient flow dynamics of dislocations which results in a fractional reaction-diffusion equation (Monneau and Patrizi, 2012).

595
The well-studied jerky flow phenomenon associated with PLC band formation is now understood as the result of initially independent, statistically distributed meso-scale dislocation mills ultimately coalescing and propagating into a macroscopic band. The best analysed mechanism for these mesoscopic fluctuations is dynamic strain ageing which gives rise to an additional characteristic retardation time scale to dislocation glide. Dynamic strain ageing operates by diffusing atoms (solute clouds) that pin dislocations and temporarily arrest the glide dislocation segments. Therefore, the plastic strain rate cannot respond 600 instantaneously to changes of the stress, and a meso-scale diffusion length/time scale emerges which stabilises the slip. The key to the oscillatory behaviour is the opposite strain-rate softening effect where the effect of retardation by diffusion is overcome through the increased mobilisation of additional dislocations (increasing the reaction rate) aided by thermal activation, hence increasing the disorder and thereby the average flow stress (Zaiser and Hähner, 1997). A similar effect of competing reactiondiffusion processes has been found in many other materials such as metal alloys (Brechet and Estrin, 1996) and self-oscillating 605 polymer materials (Masuda et al., 2016). In an elegant discourse about possible earthquake nucleation mechanisms, Orowan (1960) postulated that earthquakes may indeed be triggered by the equivalent effect of an oscillatory response of creeping rocks at depth in the crust or mantle. Orowan (1960) drew the analogy of creep failure of an annealed steel where the stress oscillates between upper and lower yield stress forming oscillatory bands called "Lüders" bands. Plastic materials, in general, are capable of arresting such small-scale fluctuations when slip tends to localise the meso-scale flow localisation and strain hardening 610 ensures that the instability mechanism cannot go catastrophic. Orowan (1960) appealed therefore to the synchronising effect of thermal feedback at depths in the Earth's mantle where self-acceleration of creep is conceivable through an avalanche-like increase of deformation where shear heating occurs faster than conduction of heat from the shear plane, finally resulting in a localised melting event as an earthquake instability.
quakes such as the 1994 great Bolivian earthquake (Kanamori et al., 1998), the mechanism may be considered an end member. Kanamori and Brodsky (2001) proposed that there are a great number of other possible earthquakes micro-instabilities without invoking a melting instability, however, the critical aspect remains to identify the physics of connecting small and large scales.
If we consider THMC feedback as a source mechanism for earthquakes, we conclude that the largest scale coupling effect is again the temperature. This becomes obvious from inspecting Table 1 where the thermal diffusion process defines the largest 620 diffusional length scale. Thermal coupling is so efficient that a material can be considered macroscopically homogeneous and at thermostatic equilibrium, while at the mesoscopic scale it still shows widespread thermal fluctuations. If we refer again to the PLC-effect as an analogy of a thermally activated material it has been recognised that these statistically uncorrelated fluctuations may, for certain critical conditions, be coordinated in the shape of a thermal wavefront that propagates through the material (Zaiser and Hähner, 1997). These wavefronts are typical acceleration waves illustrated in Fig 4 of  Lieb et al., 2020). While metals and rocks at depth are thermally activated materials, we need to discuss acceleration waves in brittle materials where the thermal activation process is less obvious.

Compression of rocks and rock analogues in the laboratory
Crushable granular materials show a similar strain-rate weakening behaviour as metals. This effect may lead to the phenomenon of propagating compaction waves in crushed snow (Barraclough et al., 2017) and puffed rice (Guillard et al., 2015) where ex-630 cellent experiments are available. The generic experimental configuration is shown in Fig. 10. Travelling compaction waves that reflect from the boundaries have been recorded by optical means (Guillard et al., 2015;Barraclough et al., 2017). Experiments have also been performed by partially soaking the puffed rice at the bottom of the experimental setup (Einav and Guillard, 2018). The interference of the compaction in unsaturated compaction and capillary-driven crushing of the puffed rice led to oscillatory catastrophic events of global compaction with acoustic-emissions perceptive as loud audible beats termed 635 "rice-quakes" by the authors (Einav and Guillard, 2018).
An important aspect of the interpretation of the laboratory experiments is that, when going from a perspective of propagating chemical or thermal acceleration wavefronts to a recording of the mechanical response of the medium, a coarse-graining step is made. We are, therefore, in most cases, not able to record the expected individual multiscale THMC wavefronts but are more likely to record the cumulative, convolved effect of the waveforms that underpin and support the mechanical deformation that 640 can ultimately lead to macroscopic failure. From a mechanical perspective, we encounter two different dynamical regimes: in the first regime the dynamics of chemical, hydraulic, elastodynamic mechanical, and thermal wavefronts is still important, and the system cannot be simplified by a thermostatic assumption. The other regime is where the thermostatic assumption holds, and the problem can be described by a quasi-static framework, where the time dependence of the system can be neglected and the problem reduces to an ideal plastic time-independent one, only controlled by the kinematic boundary conditions, e.g. the 645 position of the porous ceramic platen applying the compression in Fig. 10.
The thermomechanical theory of internal variables (Maugin and Muschik, 1999) is a variance of classical plasticity theory as it allows for additional time-dependent processes describing the evolution of local thermostatic equilibrium states through internal state variables. The introduction of additional time-dependent processes leads to a relative difference from the ideal plastic kinematic reference frame, and the localisation band can move with respect to the static, ideal plastic solution which 650 is fixed concerning the kinematic boundary conditions in a self-similar way. We can identify the two theoretical end-member cases by recording the dependence of the velocity of the optically recorded compaction waves on the velocity applied to the boundary of the experiments. If there is a positive correlation with the velocity applied to the boundary, the thermostatic endmember applies. Whereas if the velocity of the wave proves to be a material constant and independent of the boundary velocity, the dynamic solution applies. This property can be used as a diagnostic tool for determining dissipative material properties as 655 described in the appendix.

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(Dieterich, , 1987Ruina, 1983;Tse and Rice, 1986;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 diffusion matrix

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). Cross-diffusion, in excitable media (Molotkov and Vakulenko, 1993), can lead to the formation of self-supporting standing wave solutions (Berenstein and Beta, 2012). The integrable focusing non-linear 670 Schrödinger equation bears similarities but cross-diffusion waves describe a much broader class of waves and hence offer a richer solution for real-life applications (Tsyganov and Biktashev, 2014). It is, however, difficult to derive exact analytical propagating wave solutions of the cross-diffusion wave phenomenon. The piecewise linearised FitzHugh-Nagumo oscillator (Zemskov et al., 2017) is by far the best analysed prototype of a cross-diffusion wave triggered by an excitable system with short lived-spikes and a slow recovery. The FitzHugh-Nagumo cross-diffusion model has originally been developed to describe 675 the propagation of electrical impulses in nerves (Antonioletti et al., 2017).
We use a generic application of this oscillator to illustrate the potential feedback of just two THMC processes and the conditions that may prepare the formation of a future fault plane. The real waveforms of the transient cross-diffusion waves are expected to be more complex, and further work is required for a full theoretical assessment of cross-diffusion waves (Tsyganov and Biktashev, 2014). We, therefore, present only a high-level application of wave solutions to transfer knowledge 680 from other disciplines to this new field of research (e.g. mathematical biology, computational chemistry, ocean and ice waves, and photonics). We use analytical solutions developed in these fields, supplemented by our long-wavelength solution for the hydromechanically coupled case (Veveakis and Regenauer-Lieb, 2015). This case is currently being analysed in detailed laboratory experiments for testing the prediction of the cross-diffusion wave hypothesis. We, therefore, restrict ourselves here only to schematic illustrations and focus on an in-depth discussion of how the cross-diffusion concept may be applied to the 685 earthquake phenomenon.

Cross-diffusion waveforms
THMC cross-diffusion waves discussed in this work stem from multiscale fluctuations as possible wave sources with the superposition of multiple waves in a wide frequency spectrum. The convolution of these waves results in interesting dispersion patterns that bear similar characteristics to quasi-solitons encountered in optical systems where chromatic dispersion is strong 690 leading to anomalous dispersion patterns that, unlike solitons, come in discrete portions (Paschotta, 2008). Quasi-solitons encountered in photonics have two oscillating tails, one going to the left with a different wavenumber to the one going to the right. If the amplitudes of the tails are small, quasi-solitons can be treated as slowly decaying real solitons which lose their energy by radiation to form the tail (Zakharov et al., 2004).
recently (Tsyganov and Biktashev, 2014), and a broader definition of the quasi-soliton wave has been adopted to cover the richer field of wave interactions encountered in nature. Quasi-solitons feature complex dispersion relationships, where the wave velocity of individual waves (phase velocity) have a different velocity to their smooth envelope wave groups (group velocity). This leads to dominant wave packet solutions akin to envelope solitons in the nonlinear Schrödinger equation shown to be similar to the linear cross-diffusion equation in equation (31) (Tsyganov and Biktashev, 2014). They differ, however, 700 from the classical soliton solution in that the shape and speed of such waves in the established regime do not depend on initial and boundary conditions and are 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 with an approximately uniform group velocity. 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 705 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 (including 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 their original form gradually. This 710 is another feature that is different from true solitons which do not change on impact. In two-dimensional systems, additional complexities arise as they can penetrate, break apart on collision, or reflect into different directions (Tsyganov and Biktashev, 2014). Zakharov et al. (2004) compare quasi-solitons with unstable particles in nuclear physics. This describes their behaviour on collision but does not capture their capability to coordinate into standing waves on long time scales (Regenauer-Lieb et al., 2020). We propose that constructive collision of these waves may, in geological applications, lead to the bar code imprinted in 715 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., 2020) for cases where the wave collides and collapse after the 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 acoustic emission when annihilating after collision with the top surface. We propose that for a homogeneous plastic zone as shown 720 in Figure (??), the first wave-collisions on the symmetry axis of the future fault plane converts the energy loss of the crossdiffusion 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.
Another important aspect is the cross-scale coupling of the full spectrum of diffusion waves which may be seen as an energy 725 cascade from extremely short wavelength chemical cross-diffusion waves to the longest wavelength which in most cases is modulated by the thermal and mechanical diffusion length scale as shown in Table 1. The spectral content of THMC wave interactions is extremely rich, and we have emphasized that it is perhaps best compared to a convolution sharpening filter for cases of instabilities where wave energy focusses on specific locations (Regenauer-Lieb et al., 2020).
The wavelength of chemical cross-diffusion waves is relatively short and because of their rapid decay, the cross-diffusion waves emanating from either side of the elastoplastic boundary may never collide in the centre. The largest wave energy that can be transferred by cross-diffusion waves into the future fault plane is expected when the entire THMC sequence cross-diffusion phenomenon is triggered over eight to nine orders of magnitude for serpentinite ETS instability which can be modelled by the Lotka-Volterra type model (Poulet et al., 2014b).

Cross-diffusion collisions
In the context of an earthquake instability, shorter wavelength cross-diffusional waves might cascade to a small dispersion limit of a coordinated long-wavelength instability around a dominant wavenumber with a maximum amplification factor R(k) (see Eq. 26). Recently, a cascade mechanism for coupling seven orders of the length scale of quasi-solitary waves has been discovered in optical fibres (Eberhard et al., 2017a). This resonant-like scattering mode may arise through incident crossdiffusional waves colliding with reflected cross-diffusional waves in a finite width around the future fault core. On the grounds 740 of the above theoretical considerations, we postulate a series of evolutionary steps that may lead to a macro-scale rogue-wave earthquake instability. These are: -Step 1 (illustrated in Figure 11); 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 for cross-diffusion waves. In terms of anology to earthquake events the plot is understood as a convolution of the uid and solid pressure shown in Step 2 into a single curve of wave height. Note we also switch from a space representation to a time-plot on the fault plane. Step 1 (t = 0) Step 2 (critical condition for waves met) Step 3 (earthquake event, on fault plane) Step 4 (possible postseismic pattern) Wave height record of the Draupner Platform Rogue Wave event Normalized width of fault damage zone Establishing a plastic zone for the future fault damage zone Direction of traveling cross-diffusion waves from both sides Figure 11. Proposed earthquake cycle: Step 1) The size of the plastic zone is dependent on the boundary conditions and the micro-physics accommodating plastic deformation. We assume in our discussion that the zone is deforming pervasively by a power-law.
Step 2) For critical conditions (R(k) > 0) excitable self-oscillatory cross-diffusion waves propagate from opposite sides towards the future fault plane. Here, we use the analytical solution of the piecewise linear approximation of the FitzHugh-Nagumo cross-diffusion wave (Zemskov et al., 2017) for illustration. A snapshot is shown just before the collision on the future fault plane. The pressure wave of the two phases C1 (fluid, black) and C2 (solid, blue) is shown in the graph. The collision of the waves can lead to an extreme sharpening of the waveform collecting wave-power over eight to nine orders of magnitude into a rogue wave event.
Step 3) Rogue wave event on the Draupner Platform in the North Sea (Cavaleri et al., 2016).
Step 4) The cross-diffusion waves can turn into real solitons if the amplitude of their radiative tails becomes vanishingly small. We have introduced a new approach to THMC instabilities using the concept of cross-diffusion. To simplify the application, we have so far only focussed on dissipative cross-diffusional P -waves. We may now discuss the inclusion of cross-diffusional S -waves and speculate about the interaction with the P -waves. We have shown (Regenauer-Lieb et al., 2020) that the wave equations can be decomposed by using the Helmholtz decomposition of the equation of motion, into P − (compressional) and S− 765 (shear) cross-diffusion waves. 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, a bi-directional energy cascade of acceleration cross-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) 770 shear-wave instability.
The multiscale THMC-waves initially are low-energy release volumetric diffusion P -waves shown in Step 2 (Fig. 11) and are free of kinetic energy, but if 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 775 counterparts (Cartwright et al., 1997), and an earthquake instability occurs. The rogue cross-diffusion P -wave may thus be understood as a trigger to lubricate the earthquake fault. The cross-diffusion S-wave may by itself develop a similar cascade of energy and create its rogue wave, which is the earthquake instability itself.
We have also presented in Step 3 (Fig. 11) the iconic Draupner Wave example recorded on 1 January 1995 by a downward looking laser beam on the oil platform. The exceptional nature of these rogue waves bears many similarities with earthquake 780 events in terms of predictability. A hindcast analysis of the Draupner wave came to the conclusion that rogue waves are not predictable with current statistical means (Cavaleri et al., 2016). Records at a specific location are often misleading and the probability of detecting a rogue wave must be considered both in space and time; they have thought to be accepted as part of the reality of the ocean at a given sea state (Cavaleri et al., 2016). Recent work refutes this conclusion and raises hopes that a statistical analysis is possible based on large deviation theory combined with the simplified non-linear Schröedinger equation 785 solution of the cross-diffusion waves using random initial data (Dematteis et al., 2019). This raises hopes for earthquake research, acknowledging that it took nearly 15 years of research on the Draupner Wave to come to this point. To make progress on this matter, we have to leave the analytical assessment of the problem presented in this paper aside and investigate numerical solutions.
Numerical solutions of this stiff problem are unfortunately difficult, and it may be useful to consider weak-formulations 790 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). Step 4 (Fig. 11) can be obtained (Wang et al., 2018).
The quasistatic spatially inhomogeneous solutions of Step 4 (Fig. 11) are also known as Turing patterns (Turing, 1952). The special role of cross diffusion for the emergence of Turing patterns is their capability to generate spatial periodicity (Berenstein 800 and Beta, 2012), an aspect that is often overlooked in the analysis of Turing patterns. The upscaled Korteweg-de Vries equation does not necessarily honour the co-dependence between the diffusion and reaction rate constants of the cross-diffusion process, and results need to be investigated with caution (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 grow exponentially with time. To regularize the problem, finely tuned non-linear terms need to be introduced in the reaction-diffusion equation

805
(such as in the nonlinear Schrödinger equation) or, more conveniently, linear cross-diffusion-like terms need to be added in a coupled system of equations. 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.

810
Further theoretical considerations on the nature of cross-diffusion waves may also 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 an impact while the dynamics of quasi-solitons on impact are 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 815 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.
Appendices 820 A Derivation of dissipative material properties from laboratory measurements The velocity measurement in the laboratory requires an identification of a relative reference frame. When shifting from a material reference frame inside the rock formation, as discussed before, to a laboratory or field reference frame we have to consider the relative direction of the travelling wave and subtract or add the velocity vector of the boundaries of the experiment to calculate the material velocity vector of the wave from a video recording. In practice, the material wave velocity is expected to be very much faster than the background geological strain rate in the field or the velocity of the boundary, and the correction is only a minor one and may be neglected (Guillard et al., 2015). However, there exists a theoretical limit of vanishing wave speed of acceleration waves in which this is important (Barraclough et al., 2017). This limit underpins the classic theory of localisation phenomena in plasticity (Rudnicki and Rice, 1975), where the acceleration wave speed is set to zero (standing soliton wave) and acceleration waves transform into stationary localisation bands. Although, as emphasized before, the comparison with 830 conservative systems is inappropriate, we can compare this situation with a resonance phenomenon in elasticity. In the case of a diffusion wave, this situation is equivalent to the diffusion wave having the capability to retain the quickly fading information of the boundary conditions and convolve constructively with the reflected wave from the opposite boundary. This standing soliton wave solution is only possible for long wavelengths, low-frequency end members due to the physics of high-frequency damping. The standing wave phenomenon is identified as the final case of a (thermostatic) travelling wave solution with zero 835 velocity.
If localisation phenomena in the form of travelling waves or standing waves can be observed, a new diagnostic test for deriving dissipative material properties from laboratory experiments can be suggested. This test reveals how close the experiment is to thermostatic equilibrium. If the wave is controlled by its own dynamics and no relationship between the applied background strain rate of the experiment (e.g. velocity of the piston) applies, the wave velocity can be approximated by the square root 840 relationship of the Fisher-Kolmogorov type (Eq. 13). The reason for the independence of dynamic waves from the boundary velocity is that the waves quickly loose the information of the initial conditions imposed by the step function applied at the boundaries of the experiment. These travelling wave solution are a new type of waves having an extremely rich dynamics. They are quasi-solitons with very different properties to the classical soliton waves (Tsyganov and Biktashev, 2014). On the other hand, for the case of near thermostatic equilibrium, this relationship should change to a dependence on the applied bound-845 ary velocity, as the information of the applied boundary condition can spread through the sample and the diffusion wave has achieved a constructive convolution of the waves. An extreme quasi-static, thermostatic solution is shown in Step 4 (Fig. 11) where a stationary (standing) wave is shown with zero velocity. The derivation of dissipative material parameters for the zero velocity case has already been discussed (Regenauer-Lieb et al., 2013a). For the thermostatic case the wave velocity should be zero or be dependent on the boundary velocity, as discussed next. This suggested diagnostic tool opens a new way to quantify 850 dynamic material properties and to identify how close the system is to thermostatic equilibrium.

A.1 Thermostatic soliton waves
Classical plasticity theory can be extended to account for the slow dynamics before approaching the full quasi-static case.
For this case, we may neglect the dynamics of forward and backward reflecting waves from the boundaries and introduce an internal state variable that captures the relaxation to quasi-static equilibrium (Jacquey and Regenauer-Lieb, 2020). The 855 dynamics are assessed by assuming a meso-scale thermostatic equilibrium state that is linked to the next local equilibrium state by a finite-rate process slaved to the time scale of the irreversible processes (Maugin and Muschik, 1999). This is captured by the characteristic dissipative time scale of the selected internal variable α which defines the relaxation to the equilibrium of the time-dependent processes. This time scale may be defined as 860 and the velocity of the thermodynamic wave of internal variables may hence be identified as For the thermostatic case the experiment records a wave velocity that is dependent on the applied boundary velocity and we can use the above discussed theory of internal variables. For this case material parameters can be derived through measuring the diffusive length scale L d (e.g. the thickness of the propagating compaction band) and the background strain rate asα (e.g.

865
compaction strain-rate imposed on the sample by the compression piston). This interpretation assumes that all deformation activity is taking place in the active band. We can invert the magnitude of the internal variable α (e.g. the compaction strain) from measuring the material speed of the propagating compaction band by using the mass conservation criterion in Eq. 35. This method of deriving the band strain in a propagating compaction band in compaction of snow has been used by Barraclough et al. (2017). 870 We have performed a similar compaction experiment with a highly porous limestone which produces first a stationary compaction band roughly in the middle of the sample after only 3% axial strain (Chen et al., 2020). When interpreting the result in terms of a standing wave limit (zero acceleration wave speed) we expected the location to be strictly linked to the quasistatic geometric constraints which can be either the size of the sample or an internal defect. Our time-lapse X-ray Microscopy experiment (see Chen et al. (2020) for details) revealed a local porosity variation as the initial seed of the compaction band.

875
In subsequent stages of the experiment, the stationary discrete compaction band in turn acted as the nucleation point for a compaction wave travelling downwards from the discrete band ultimately leading to full compaction of the lower part of the sample. This observation is strictly identical to the experiments performed by Barraclough et al. (2017).

A.2 Dynamic quasi-soliton waves
In the case where the experiment records a wave velocity that is independent of the applied boundary velocity, we can derive the 880 dissipative material properties directly from measuring the diffusion length scale L d and with the measured material velocity inverting Eq. 13 to obtain the effective Damköhler number Da. Alternatively, we can derive the instantaneous modulus C from inverting Eq. 32 from a measurement of the density ρ and the band velocity v. Guillard et al. (2015) have performed experiments with compacted puffed rice in a setup similar to Fig. 10 recording the alternative case of a dynamic solution. The wave speed was found to be independent of the applied boundary velocity. Guillard 885 et al. (2015) found that, when varying the velocity by a factor of 2.5, the number of propagating bands increased but not their velocity. If the density of the travelling compaction band can be measured by high-speed (Synchrotron) X-Ray-CT or radiography the instantaneous modulus C can directly be inverted from the band velocity. Such a measurement could then be used to derive the effective elastodynamic modulus quantifying the visco-elasto-plastic properties of the material. Kolmogorov, A., Petrovsky, I., and Piskunov, N.: Etude de l'equation de la diffusion avec croissance de la quantite de matiere et son applica-