Cross-Diffusion Waves as a trigger for Multiscale, Multi-Physics Instabilities: Theory

. We propose a non-local, meso-scale approach for coupling multiphysics processes across the scales. The physics is based on discrete phenomena, triggered by local Thermo-Hydro-Mechano-Chemical (THMC) instabilities, that cause cross-diffusion (quasi-soliton) acceleration waves. These waves nucleate when the overall stress ﬁeld is incompatible with accelerations from local feedbacks of generalized THMC thermodynamic forces that trigger generalized thermodynamic ﬂuxes of another kind. Cross-diffusion terms in the 4 × 4 THMC diffusion matrix are shown to lead to multiple diffusional P - and 5 S -wave equations as coupled THMC solutions. Uncertainties in the location of meso-scale material instabilities are captured by a wave-scale correlation of probability amplitudes. Cross-diffusional waves have unusual dispersion patterns and, although they assume a solitary state, do not behave like solitons but show complex interactions when they collide. Their characteristic wavenumber and constant speed deﬁne mesoscopic internal material time-space relations entirely deﬁned by the coefﬁcients of the coupled THMC reaction-cross-diffusion equations. Part 2 proposes an application to earthquakes showing that for ex-10 treme conditions, cross-diffusion waves can lead to an energy cascade connecting large and small-scales and cause solid-state turbulence


Introduction
The theory presented in this paper grew out of the conference series dedicated to understanding Coupled Thermo-Hydro-Mechanical-Chemical (THMC) in Geosystems (GEOPROC).The 7th international event was held in 2019 in Utrecht and focussed on earthquake and faulting mechanics (this special volume).Integration of mechanical, hydrodynamical, thermal, and chemical processes covers, however, a much wider field from the pore to plate tectonic scale for a wide range of natural and engineering problems in geological system discussed in focus topics on earlier GEOPROC conferences.These problems include nuclear waste disposal, coal seam gas, enhanced oil and gas recovery, geothermal energy, mineral deposits, tailing dam collapse, landslide and many others.The individual problems may have their own characteristics.However, the common scientific issue of multi-scale feedback of THMC processes remains the same.
The GEOPROC theme seeks to foster the urgently needed growth of experimental, numerical, and theoretical studies on multiphysics (THMC) and multiscale framework studies in Earth Sciences.The current practice is to still use engineering solutions based on empirical material laws to address specific natural and engineering problems in geological systems and energy production in geothermal energy, nuclear waste disposal, reservoir engineering for oil and gas, the formation of mineral deposits, induced seismicity, natural hazards, and CO 2 sequestration and utilization.These empirical engineering approaches are often inadequate, as indicated, for example, in the failure to avoid the 5.5 magnitude earthquake in Pohang Korea in November 2017, which was anthropogenically induced by high-pressure hydraulic injection during the previous two years (Grigoli et al., 2018).
Part of the reasons for lack of a wider adoption of coupled THMC approaches in the community is a lack of a theoretical basis on which to assess the rich solution space that arise from a coupling of the four (THMC) partial differential reaction diffusion equations.While parallel numerical tools for modelling fully coupled non-linear systems of THMC equations have become available through pioneering work in nuclear engineering (Gaston et al., 2009;Permann et al., 2020), the corresponding theory has not progressed as far.The application of the powerful nuclear engineering modelling tool has been successfully transferred to geosciences and applied to geodynamic modelling (Jacquey and Cacace, 2020a, b) and the modelling of the Non Volcanic Tremor and Slip (NVTS) events in the circum-Pacific subduction zones (Poulet et al., 2014b) as well as applied to geological faulting problems (Poulet et al., 2014a).However, a sound theoretical description and interpretation of the local processes resulting in the interesting macroscopic phenomena has been lacking.Part 2 (Regenauer-Lieb et al., 2020) will attempt to rectify this shortcoming.
Before discussing a possible application of the new theory to the processes of earthquakes and faulting in our companion paper (Regenauer-Lieb et al., 2020), here we present a transdisciplinary approach bridging the gap between observations of instabilities from the molecular scale to the very large scale.The theory in this paper is written using approaches familiar to the theoretical and applied mechanics community.The original work is based on the 1960's work (Hill, 1962) building the foundation of theoretical approaches to localisation criteria, via the so-called acoustic tensor criterion, widely used in the engineering community (Rudnicki and Rice, 1975).The approach focusses on standing wave quasi-static solutions based on vanishing speeds of acceleration waves.Surprisingly, little effort has been made to explore the rich wave field of the corresponding travelling wave solutions, probably because dynamic events are only of academic interest to the engineering plasticity community that focusses mainly on developing safety standards as well as limit analysis and design.A notable exception is the work of Benallal and Bigoni (2004) who found that under dynamic conditions, unbounded growth of perturbations can be found in the short wavelength regime with divergence growth.This calls for an extension to the theoretical work of Hill (1962) which is presented here.
The dynamic field is, however, of special interest to the researcher in the area of earthquake and faulting instabilities.The state of the art in this field is defined by the influential experimental work of Dieterich (1979) including the early work on the application of the rate and state variable friction approach to earthquakes (Tse and Rice, 1986).The approach based on these laboratory-derived constitutive equations has reached a mature stage, and no attempt is made here to compare the rich field of findings with the present theory.We approach the problem from an entirely different angle through theoretical investigation of the mathematical solutions of the system of coupled partial differential THMC equations that deliver wave solutions with short wavelength instabilities.In the course of developing the new approach, we describe wave physics phenomena that have previously not been reported in geophysical solids but are well known in a range of different fields from quantum systems to ocean waves (Zakharov et al., 2004).It is fair to say that the theory is rather in its infancy state, and special care needs to be taken before considering a direct application to the aforementioned systems.The first part therefore presents the theoretical derivation, and the second part delves into possible applications and proposed experimental tests to verify the applicability of the theory.
In this paper, we introduce the classical approach of acceleration waves in plasticity theory to the seismology community by starting with the Helmholtz decomposition of the seismic wave equation into P-and S-waves (see section 3.1).We show how plasticity theory can be integrated via Hill's acceleration waves into the equations.This approach leads directly to the unbounded short wavelength growth described by Benallal and Bigoni (2004) which cannot be solved without further assumptions.The innovation proposed in the two manuscripts is to appeal to the multi-scale nature of the THMC coupled problem.
We regularise the problem by embedding an open system thermodynamic meso-scale theory with unbounded solutions into a closed system macro-scale approach that describes the emergence of a standing-wave solutions.
There are two opposite starting points for the derivation of the approach.Here, we investigate the meso-scale from the conventional mechanical quasi-steady state (infinite time scale) solution of the macro-scale.The present manuscript uses the macro-scale view, i.e. the classical mechanical viewpoint for the investigation of the physics of acceleration waves (Hill, 1962).
The second part (Regenauer-Lieb et al., 2020) describes the meso-scale view which is the classical viewpoint of a chemist.
Both viewpoints are objective descriptions of the coupled THMC problem and should deliver the same outcome.
In order to define the separation between the meso-and macro-scale of a THMC coupled problem, we propose that the scale for each of the THMC-processes is defined by their own characteristic diffusion time/length scales (Regenauer-Lieb et al., 2013b).For simple problems progress can be made by studying thermodynamic equilibrium states in isolated closed systems.Likewise closed, coupled, far-from-equilibrium THMC systemss that feature irreversible behaviours can be modelled by a thermomechanics approach (Collins and Houlsby, 1997), also called a thermodynamics with internal variables approach (Maugin and Muschik, 1999) or a hyperplastic approach (Houlsby and Puzrin, 2007).This theory is, however, only applicable to faults that have reached a thermal steady state as implied by a standing wave solution of acceleration waves.This approach prevents modelling of dynamic phenomena.Modelling of earthquake and faulting is hence one of the most difficult topics to address using a self-consistent thermodynamic approach.
A particular challenge for deriving dynamic THMC coupled wave solutions is the discrete nature of the cascade of steady state solutions defined by the standing wave solutions of thermomechanics which leads to a discrete material behaviour as discussed in the next section.Standard probability theory is therefore not suitable as this assumes a continuum of wave functions (Cohen, 1988).In order to solve this issue we use a transfer of knowledge from classical quantum mechanics to characterise any system at a larger scale.The information on multiple internal material time/length scale processes travels and disperse each at characteristic velocities in the form of acceleration waves.
The nucleation mechanism of these waves relies on the meso-scale open system behaviour where the overall macro-scale thermodynamic forces can become incompatible with accelerations from local thermodynamic fluxes.These incompatibilities radiate wave energy away from its source in the form of "cross-diffusion" waves.The emergence of cross-diffusion waves can be perhaps best understood from a chemical viewpoint (Regenauer-Lieb et al., 2020) where propagating chemical waves have been studied in details (Vanag and Epstein, 2009).In chemical systems cross-diffusion, is defined as the phenomenon in which a gradient in the concentration of one species induces a flux of another chemical species.In the present context thermodynamic forces and fluxes are generalized THMC fluxes defined in Table 1.Before discussing cross-scale coupling of thermodynamic forces and fluxes in sections (3.3) following, it is useful to briefly review insights into the formation of discrete dissipative structures.

Dissipative structures
The concept was introduced first in chemical and biological systems where morphogenic patterns (Turing, 1952) were identified as solutions to the underlying reaction-diffusion equations (see Figure 1).These discrete patterns were later on named Turing patterns.A review of Turing patterns in nature can be found in Ball (2012).
We propose here a generalised approach to cross-diffusion that is known in bio-physics as taxis (Heilmann et al., 2018).For example, the pufferfish and the siltstone in Figure 1 show similar patterns caused by a fundamental mechanism known as taxis.This is a process that forces components of a pattern to organise as an ensemble in reaction to changes in the environment, and so to move towards, or away from a perturbation .The process eventually leads to the formation of a new energetically stable pattern.The fundamental pattern-forming taxis mechanism can be caused by adhesive forces (hapto-taxis), hydrodynamic (gyro-taxis), gravitational (gravi-taxis), light intensities (photo-taxis), or chemical driving forces; chemo-taxis; as shown in Figure 1.Non-biological patterns are generally formed by Thermal, Hydrodynamic, Mechanical, and Chemical (THMC) reactions, which can also include electrical and biological drivers.In mathematical biology, taxis models are used to understand and quantify a variety of complex problems, ranging from nerve pulse responses and spreading of diseases (Zemskov et al., 2017) to predicting the spatio-temporal patterns of predator-prey systems.
However, except for the seminal early work by Ortoleva and co-workers (Dewers and Ortoleva, 1990;Ortoleva, 1993Ortoleva, , 1994) ) developments of taxis models in Earth and Material Sciences have lagged.A review of the progress made in this field as well as a specific case study of rhythmic banding in marls can be found in the recent work of L'Heureux (L'Heureux, 2018(L'Heureux, , 2013)).
By analogy to the mathematically similar biological and chemical systems we propose here that earthquake instabilities are preceded by propagating THMC dissipative waves which could enable new detection methods if they can be resolved by sensors.We will discuss such possible precursor phenomena for earthquakes in part 2 (Regenauer-Lieb et al., 2020).In chemical systems, propagating waves stemming from reaction-diffusion physics, are very well documented and the hypothesis Both pattern formations can be modelled by a reaction-diffusion-advection type instability due to chemo-taxis.
here is that this applies generally to all THMC coupled processes.Excellent update of the formulation for chemical systems can be found in Vanag and Epstein (2009).While chemical oscillations thus appear to be well understood, the phenomenon of an oscillatory response is less well established in other THMC systems.However, these systems show the same transitions from a simple continuum response to a highly localized state.The existence of a discrete, particle-like nature has also been discovered in fluids when they are driven far from equilibrium.If driven far from equilibrium by surface forces, fluids clearly show (see Fig. 2) a highly dissipative, sharp transition from a continuum state to one of a highly localized, propagating, particle-like state (Lioubashevski et al., 1996).
For the case of deforming geomaterials localization phenomena, characterized by a sudden transition from continuum deformation behaviour to a highly localized state, are also well established in solid mechanics applications (Rudnicki and Rice, 1975).Fig. 3 shows a periodic set of localised deformation structures formed as a result of a compressive tectonic regime.Similar standing wave like features are encountered in many geological systems L'Heureux (2018); Ball (2012).However, direct experimental evidence for precursory transient travelling solitary states is largely unknown and has only been shown recently based on mathematical considerations (Hu et al., 2020).In the supplementary material of part 2 (Regenauer-Lieb et al., 2020) we will discuss possible experimental tests of the precursor phenomena.The lack of experimental evidence can be explained by the challenging task of dealing with the large length scale of the geomechanical phenomena and the long-time scales of observation required to mimic natural processes in the laboratory.
The dynamics of the formation of these mechanical dissipative patterns can therefore only be investigated using analogue materials in the laboratory.Analogue experiments have been performed in a granular, brittle-matter compressed in a uniaxial direction (Guillard et al., 2015;Einav and Guillard, 2018).A propagating compaction wave phenomenon has been observed.into acoustic emissions.The phenomenon has been compared to ice-quakes in ice sheets (Einav and Guillard, 2018).While these experiments allow some insight into the precursor phenomena of stationary compaction bands the experiments themselves never reached the stationary mode.This aspect will be discussed in more detail in part 2 (Regenauer-Lieb et al., 2020) where experiments that are very close to the stationary mode are also introduced (Barraclough et al., 2017).The stationary mode allows development of a robust thermomechanics theory which offers a modular thermodynamically self-consistent approach for modeling earth instabilities (Jacquey and Regenauer-Lieb, 2020).

Acoustic bursts have been registered when the waves interact with interfaces leading to the conversion of their kinetic energy
Experiments with highly porous carbonates have been performed (Chen et al., 2020).These produced stationary and nonstationary compaction bands under uniaxial loading.Unfortunately, these experiments have the opposite problem in that the exact analysis of the dynamic evolution did not reach sufficient resolution in space and time to convincingly detect the wave phenomenon described in the granular brittle matter.We, therefore, explore in this contribution theoretical predictions of the dynamic wave propagation based on an extension of the thermomechanics approach for wave propagation in dissipative materials (Coleman et al., 1965).
We propose here that the dissipative wave phenomenon is universal for THMC reaction-diffusion systems that are driven far from equilibrium.The approach allows an interpretation of observations in nature and the laboratory in terms of propagating particle-like states which emerge as stationary Turing patterns for long-timescale standing wave solution of a THMC cross- diffusion formulation.In order to recover the dissipative wave equations, we present in the following the standard constitutive assumptions for any generic thermodynamic fluid or solid mechanical system and describe how the physics of THMC feedbacks can be implemented to resolve the phenomenon of propagating dissipative waves in these systems.

Constitutive assumptions
The fundamental equation of motion is: where σ = σ ij is the Cauchy stress tensor, ρ the density, f is a body force (e.g.gravity) and a the acceleration.This equation does not stipulate a constitutive law but with constitutive assumptions it becomes the master equation for the theory of elastic waves, fluid mechanics and continuum mechanics.In elasticity, wave equations directly result from the equation of motion defining the wave characteristics by using the Helmholtz decomposition, in terms of shear (S-wave) and compressional (Pwave) wave velocities which is a convenient description for the purpose of this paper.
For an isotropic elastic medium, for instance, accelerations in Eq. ( 1) are only allowing elastic displacements described by u.In this case the material can be characterized by just two velocities: the elastic P -wave velocity v p and the elastic S-wave velocity v s , and we obtain from Eq. ( 1) the elastic-wave equation: . (2) Similarly, by allowing the material to deform in a viscous manner, acceleration can be monitored by a local change in velocity v, and the Helmholtz decomposition identifies a scalar P -wave and a vectorial S-wave potential field.The material constants are the dynamic shear η and bulk ζ viscosities to obtain the generalized Navier-Stokes equation: where is the deviatoric viscous strain-rate.with where I is the identity matrix.
We emphasize here, that although the Helmholtz decomposition can be performed in a similar way to derive volumetric and shear moduli that describe dissipative material behaviour, their response to infinitesimal perturbation is generally to dampen propagating elastic waves.One could, therefore, come to the erroneous conclusion that their contribution to precursory wave phenomena to macroscopic failure is an overall suppression of instabilities.Coleman and Gurtin (1965) have shown that this conclusion is wrong using the concept of materials with fading memory conceptualized by rational thermodynamics.We are using a simpler approach and are introducing fading memory through THMC dissipation processes based on the non-equilibrium thermodynamics approach of deGroot (1962).Therein, non-equilibrium conditions are seen as a concatenation of thermostatic equilibrium processes.We, therefore, can use the local equilibrium definition of the pressure as p = − ∂U ∂V , where U is the internal energy and V the volume and explore the emergent dynamics through investigating the stability of small perturbations from individual equilibrium states.In the present context pressure is therefore defined as p = 1 3 tr(σ) and is negative for compression.
For the elasto-visco-plastic case, we have the equivalent fourth-order elastic-viscoplastic stiffness tensor C characterizing material stiffness and the corresponding elasto-viscoplastic bulk viscosity ζ to give in this case ˙ denotes the deviatoric strain rate which in the purely elastic case before yield is ˙ = ˙ e becoming post-yield the elasto-viscoplastic strain-rate defined by ˙ = ˙ e + ˙ vp where the subscripts e and vp refer to the elastic and viscoplastic components.The same approach is used for decomposing the equivalent elasto-viscoplastic volumetric strain-rate ˙ 0 .While the emergence of elastic P − and Swaves for any infinitesimal disturbance is a well-known physics phenomenon, dissipative processes are commonly known to dampen elastic waves.Long-range wave propagation in dissipative materials was therefore contested for a long time (Coleman and Gurtin, 1965).Ideal elastic waves without damping propagate without loss of energy as they are based on the conservation of energy.However, the dissipative chemical and biological diffusion waves are known to propagate as autonomous wave sources, spontaneous oscillations, and quasi-stochastic waves which are synchronized over the entire space to form dissipative structures (Vasil'ev, 1979).
They are an entirely different class of waves as they are based on dissipation in active kinetic systems in contrast to waves in conservative systems.For simplicity, we only discuss the slow visco-plastic wave phenomenon allowing the investigation of conservative and dissipative waves as different processes.For decoupling elastic and dissipative waves, we need to assume large differences in the propagation speed of the waves.This is done by assuming Maxwellian rheology, implying a separation of elastic and visco-plastic wave time-scales in the context of an additive strain-rate decomposition of Eq. ( 2) and (4).To recover dissipative waves from the above discussed Navier-Stokes equation, modified for the inclusion of elastoplastic behaviour, we introduce local thermodynamic THMC feedback processes that change the instantaneous fourth-order elastic-visco-plastic stiffness tensor C.

Acceleration waves and classical theories of localisation
The simplest implementation of the non-equilibrium approach of deGroot ( 1962) is the theory of internal variable thermodynamics which unifies the kinetic description of chemical reaction-diffusion processes and the above described elasticviscoplastic formulation (Maugin and Muschik, 1999).Perturbations to the local equilibrium assumption of the non-equilibrium thermodynamic theory of internal variables can lead to conditions of violation of smoothness on surfaces in a body, where one or more internal variables from the lower scale suffer jump discontinuities, owing to locally reaching a critical dissipation.In a classical thermodynamic sense, this can be viewed as suddenly switching on a micro-engine somewhere in the system that disturbs the overall stress field.This is the physical reason for the formation of acceleration fronts, where the diffusive lengthscales are linked to the convective velocity of the step function on dissipation waves.While fluid-and solid-wave phenomena occur when the above equation includes inertial forces, so-called 'acceleration waves' (Hill, 1962), caused by local surfaces of acceleration (see Fig. 4), can also occur in the creeping flow limit when no acceleration due to a (gravity) potential is present and f = 0.These acceleration waves are defined as geometric surfaces (here assumed to be plane waves) that move relative to the material.
Acceleration waves can be described in two ways.One can use two coordinate systems, one for the reference state and one for the current state.A more elegant way is to consider convective coordinate systems by formulating the constitutive law in terms of stress-rate.For this we consider the space derivative normal to the moving wavefront (see Figure 4) indicated by ∂ ∂s .Waves are travelling concerning a background Lagrangian moving material reference frame.
Considering the traction (load per unit area) in the direction normal to the wavefront as F and choosing the magnitude of the velocity of the moving wavefront as c , the jump condition indicated by the square Iverson brackets can be advected along c.This leads to Hadamard's jump condition where the true traction rate along the advected coordinates is: Hadamard's jump condition applies to all internal variables and the acceleration across the wavefront is constrained by Combining Eq. ( 5) and ( 6) we obtain Substituting v = −c ∂v ∂s into Eq.( 7) we obtain Hill's formulation of acceleration waves in Eq. ( 8) expresses the energetics of the acceleration waves by the square of the material velocity c times the mass of the characteristic segment defined by ∂v ∂s .This provides a simple formulation where the energetics of the material is solely described by Eq. ( 1) and the mesoscale mass exchange rates on acceleration waves by Eq. ( 8).The material velocity c, being the velocity of acceleration waves becomes a material constant for the propagation of acceleration waves.Acceleration waves form the basis of localization criteria in plasticity theory.The criterion for instability is derived from the equivalent theory in elastodynamics where for an elastoplastic body the acoustic tensor Γ is defined by In elastodynamics, the eigenvalues of Γ divided by the mass density represent the square of the elastic wave propagation speed in the direction of the unit normal vector n.In elasto-plasticity, the equivalent dynamic stability criterion is defined by Eq. ( 8) which in terms of acoustic tensor implies that Γ ∂v ∂s = ρc 2 ∂v ∂s (10) Dynamic system stability can be evaluated through assessing the eigenvalues of the acoustic tensor thus determining the speed of the acceleration waves which must be real and defined by the square root of the instantaneous modulus C divided by the instantaneous density ρ (Coleman and Gurtin, 1965).Mathematically, Eq. (2 -4) can be represented by the addition of two functions, a scalar field and the curl of a vector field.The former without curl or rotation identifies dissipative compressional P -waves and the latter features zero divergence and corresponds to isochoric sinistral and dextral dissipative shear S-waves.
These waves are interpreted as stationary (standing) waves when the determinant of the acoustic tensor and consequently the wave speed is zero, which is the standard condition for localisation in plasticity theory (Vardoulakis and Sulem, 1995).Accordingly, the formation of localized shear-bands out of homogeneous plastic flow is assumed, when the velocity of the wavefront vanishes.Hill (1962) was discussing shear acceleration waves in an ideal linear, time-independent elasto-plastic material where two families of characteristics (dextral and sinistral slip lines) feature a jump in strain-rate at the wavefront accompanied by one in stressrate (but not in stress).This in turn implies a related jump in stress gradient.Later work extended the theory to formulate accelerations waves as the basis of modern criteria for localization in plastic media (Rudnicki and Rice, 1975;Rice, 1976).
In those theories the possibility of volumetric acceleration waves was, however, neglected and volumetric deformation was parameterized by an empirical dilatancy angle.Another shortcoming of the localization criterion for the application to THMC instabilities is that it is not directly applicable to the rate-dependent elasto-viscoplastic case.The inclusion of rate effects implies a positive wave speed different from zero (Duszek-Perzyna and Perzyna, 1996).The thermomechanics approach (Jacquey and Regenauer-Lieb, 2020) allows incorporation of a quasi-static wave speed related to the internal variable that quantifies the rate-dependence.However, to date, no generally accepted localization criteria for the transition from a dynamic to quasistatic rate-dependent solution of Eq. ( 4) exists, although stationary and wave-like propagating localization phenomena for rate-sensitive materials (Barraclough et al., 2017) are observed in the laboratory and nature.The method of choice to date is to use all field equations and perform a numerical stability analysis.A discussion on an extension to the above-discussed criterion has been presented recently (Pisanò and Prisco, 2016) and energy-based criteria that successfully model the adiabatic limit have been repeatedly revisited many times over the past 30 years (Paesold et al., 2016).The present approach provides an alternative path to systematically analyse the full system of field equations.

THMC acceleration waves
We assume creeping flow in Eq. ( 4), and there is, therefore, no effect of gravitational acceleration (F = 0).We will show that the meso-scale formalism identifies an alternate internal force density from within the considered material volume stemming from a local thermodynamic THMC force (e.g.∇p).This internal force integrates over the accelerations a M of the microprocesses inside the continuum element by multiplying them with the average volume density.These accelerations stem from dissipative mechanisms (e.g.volume changes by phase transitions, fracture, etc.) inside the representative volume element (Hu et al., 2020).For critical conditions, they can cause acceleration waves propagating as creeping waves.Hadamard's jump conditions need to be extended for internal THMC variables µ such as temperature, porosity, permeability, viscosity, etc...
Hadamard's jump conditions state that if time derivatives ( Ḟ, v, μ) and gradients ∇F, ∇v, ∇µ) have jump discontinuities across the wavefront then F, v and µ are continuous functions of space.The compatibility condition relating jumps in rates of change of internal variables to jumps in gradients for all internal variables µ (Duszek-Perzyna and Perzyna, 1996) implying that the jump in the gradient of pressure inside the acceleration wave is constrained by Acceleration waves consider a step function (Eq.5) where the stress-rate is discontinuous along the surface.The stress is, however, continuous across the wavefront and the stress self-diffusion coefficient is also constant outside of the wave.Therefore, for modelling acceleration waves in a homogeneous material we can simplify Eq. ( 1) further and assume constant bulk and shear viscosity outside the wave and assume continuity of stress across the acceleration wave.Noting that the traction in the direction of the normal vector n on the acceleration wavefront is F = n • σ it follows from Eq. ( 7) that the jump in stress-rate on the acceleration wave is (Duszek-Perzyna and Perzyna, 1996): Substituting the stress-rate for the acceleration v from Eq. ( 13) and the pressure rate for the gradient of pressure from Eq. ( 12) and inserting the jump condition into Eq.( 4) it follows that If we define the magnitude of the wave speed in the normal reference system as w = w • n, then c = w − v • n is the local particle velocity of THMC accelerations in the acceleration wave relative to the normal material velocity.
Eq. ( 14) allows us to draw some important conclusions for acceleration waves.
(1) The first term on the right shows that the pressure rate divided by the wave velocity or the equivalent gradient of pressure plays an important role in acceleration waves.
(2) The second term on the right implies that gradients of deviatoric strain-rates are related to rate changes of the stiffness tensor as implied by the jump condition of the internal variable inside the propagating wave.Recall that the jump condition (Eq. 5 or 12) advects jumps in gradients of the internal variable around the propagating wavefront through a jump in the rate of change of the variable.
(3) The last term implies that the same is true for the volumetric strain-rates and the rate of change of bulk viscosity.

Multiscale cross-diffusion model
So far we have only discussed the mechanical reaction-diffusion equation, where the shear and bulk viscosities control the diffusion of stress.For the multiphysics implementation, it is convenient to think of diffusion of momentum and use the momentum diffusivity (kinematic viscosity) instead of the dynamic viscosity.We, therefore, denominate ζ M as the volumetric diffusion coefficient of pressure (kinematic viscosity).In the following, we first formulate the reaction-diffusion equation in a classical way.That is to say that meso-scale cross-diffusion effects are neglected.We identify THMC-Turing patterns as multiscale energy eigenstates of the reaction-diffusion equations thus characterizing Prigogine's dissipative structures if they emerge.
In these formulations, the viscous (M) mechanical pressure diffusion equation finds its counterparts in the equivalent thermal (T) Fourier-, (H) Darcy-and (C) Fick's-diffusion laws where the diffusion coefficients are indicated by the associated THMC subscript.The corresponding reaction rates are the local hidden-variable reaction rate R T , R H , R M and R C , respectively.It is common practice to ignore the meso-scale cross-diffusion kinetics introduced in the classical theories of localisation.We emphasize therefore the difference between large-scale reaction rates R i and meso-scale reaction rates r i of the non-local theory which considers the important effect of cross-diffusion.The two rates are identical in the infinite time-scale limit as cross-diffusion can be eliminated adiabatically (Biktashev and Tsyganov, 2016).
In the adiabatic limit we obtain similar reaction-diffusion equations across a vast range of THMC diffusion length scales as tabulated in ta The reaction rates most often stem from different micro-processes at lower scale inside the considered continuum element which introduces cross-scale diffusion fluxes as shown in the next section.In order to generalize the approach we propose that the composite multi-scale THMC wave operator ĤT HM C in Table 1 can be constructed through a linear superposition characterized by where ∂z 2 and i refers to the individual thermodynamic THMC process.For the following discussion, we simplify further and neglect the deviatoric terms in Eq. 14 and retain only the scalar volumetric terms and reduce the equations to 1-D.This allows us to investigate the poorly known volumetric dissipative waves which must exist in addition to the dissipative shear waves discussed by Hill (1962).To introduce the meso-scale considerations we identify wave scale reactive source terms r T , r H , r M , r C of the hidden variable rates R T , R H , R M , R C as the actual terms that trigger acceleration waves.These meso-scale source terms stem from a jump in the thermodynamic force (gradient of the variable) into a jump in the thermodynamic flux (rate of the variable, i.e. temperature, fluid, and solid pressure and concentration).The important volumetric coupling is overlooked in the classical localisation theory (Rudnicki and Rice, 1975).The wave-scale source term provides the convected pressure rate built up by internal accelerations.It relates to the local mass exchange processes according to Eq. 8.
5 Cross-diffusion as a non-local theory for localisation In the following we generalise the discussion of the meso-scale THMC mass exchange processes using mixture theory applied to HM coupling as presented in Hu et al. (2020).We show that the physics of cross-diffusion follows from a reactive source term at the macroscale that requests a cross-diffusion term at the meso-scale for thermodynamic consistency.The full derivation is found in Hu et al. (2020).Here we summarize the main conclusion from the mixture theory analysis for convenience.
We consider two mass fractions A and B for mass exchange denoted by the i th and j th phase as an example.We identify ξREV i as the large-scale Representative Elementary Volume (REV) for averaging of mass transfer rate from the phases A to B where V REV , V A , V B denote the REV volume and the volume of the i th and j th phase, respectively.ξREV i defines the REV-scale averaging of the mass exchange rate between the phases where the REV-scale source term of mass is obtained from the other species: where ξlocal Mass conservation at global scale for the phases A and B gives: ρ A and ρ B identify the density of the respective phases, while v A and v B their velocities in the direction of x while ξA and ξB represent the volume averaged mass generation in the REV.Following an approach presented in Hu et al. (2020) using the generalised THMC mass exchange processes we arrive at: In a saturated porous medium, a straightforward interpretation of ρ A and ρ B may be the density of the fluid phase and that of the solid phase, respectively (Hu et al., 2020).The interpretation of the incorporation of the effects of chemical and thermal processes may not be as straightforward for the observer as they act via Eq.15 as a linear convolution operation.If we interpret the time-domain convolution operation of THMC waves in the frequency domain, then the chemical and thermal waves can be seen as filters for HM coupling, sharpening or smoothing the waves.The interpretation of THMC waves in terms of a sharpening or smoothing filter analogue is discussed in detail in part 2 (Regenauer-Lieb et al., 2020).
To illustrate the point of choosing a particular time-space scale of observation of THMC waves we first consider the simple homo-entropic flow assumption and choose a classical mechanics point of view.Density is then defined as a function of pressure, temperature, and chemical concentrations by the Equation of State.The coupling in Eq. 19 leads to the possible nucleation of Hydro-Mechanical cross-diffusion pressure waves (Hu et al., 2020).Considering that possible thermal processes, ρ A , and ρ B , at the solid-fluid interface may be affected by the local temperature changes, we identify new time-dependent processes at the solid-fluid boundary.The process of heat transport sets two new internal timescales changing the fluid and solid pressure, respectively.Therefore, the thermal process acts as a 'convolution filter' added to the pressure evolution of each phase.Conversely, the pressure diffusion process in the solid and fluid phase triggers two additional timescales in the change of temperature.Now a 3 × 3 Thermo-Hydro-Mechanical cross-diffusion formulation can be obtained following the same steps of upscaling from local to a global scale, and the additional 4 timescales correspond to the newly introduced 4 cross-diffusion coefficients.One can arrive at a similar conclusion by considering the interplay between the evolution of pressure distribution and internal chemical reactions.

Formulation of the THMC cross-diffusion matrix
The concept of cross-diffusion is well known in chemistry.In a chemical system with just two species A and B, 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 drives another generalised thermodynamic flux.Staying with the chemical example of species A and B, we have in 1-D: where r A and r B are the local source terms using the mesoscopic self-diffusion and cross-diffusion decomposition in Eq. ( 19).
Following Eq. ( 15) we can now generalize (20) to include the full cascade of internal accelerations through multiscale coupling.
Cross-diffusion allows coupling of accelerations from one classical REV-scale reaction-diffusion system, defined by the (self- )diffusive length scale ζ (T,H,M,C) t, to another.These cross-diffusion coefficients link the gradient of a thermodynamic force C j of one THMC process to the flux of another kind.This allows the definition of a fully populated diffusion matrix as in: The cross-diffusion processes formulate the link between different THMC processes.The cross-diffusion coefficients thereby introduce new cross-scale coupling length and the time scales and are often much smaller than the self-diffusion scales.This is not always the case (Manning, 1970).Hu et al. (2020) show normal examples where cross-diffusion length/time scales are much smaller than the self-diffusion length scales.

Criterion for nucleation of cross-diffusion Waves
A detailed discussion of the criterion for nucleation of cross-diffusion waves and their waveforms can be found in Tsyganov and Biktashev (Tsyganov and Biktashev, 2014).Here, we first summarize the basic method that is well established in the fields of mathematical biology and chemistry and follow on with a discussion of other communities, where the phenomenon of cross-diffusion waves is well-documented under different names.
The criterion for nucleation of cross-diffusion waves relies on assessing the dispersion relation of the eigenvalues of the characteristic matrix of a perturbed cross-diffusion-reaction equation (Vanag and Epstein, 2009).The eigenvalues are functions of the square of the wavenumber of the perturbed state and identify the growth rate of the perturbations.This approach for deriving the mathematical criterion for nucleation of acceleration waves is hence evaluated from a small plane-wave -perturbation of Eq. ( 21) with The characteristic matrix of the thus perturbed Eq. ( 21) allows assessing the stability of the system.Accordingly, all eigenvalues of the characteristic matrix must be real and positive, and hence the determinant of the matrix must be larger than zero.For determinants smaller than zero, cross-diffusion waves are expected to propagate as quasi-solitons (Tsyganov and Biktashev, 2014).A worked example for hydromechanical cross-diffusion waves can be found in Hu et al. (2020).
6 Soliton versus quasi-soliton solutions Since cross-diffusion waves in geomaterials are largely unexplored due to the extreme length and time scales encountered in a geosystem, an appreciation of their complex characteristics can be obtained from mathematically similar systems such as waves in oceans, Lasers, and ice.There is an important difference between solitonic waves and quasi-solitonic cross-diffusion waves.We follow Zakharov et al.'s (Zakharov and Kuznetsov, 1998) definition of solitons and quasi-solitons and identify solutions to the perturbed Eq. 22 of the type: Quasi-solitons appear as non-local solutions to Eq. ( 22) if localized coherent structures defined by true soliton solutions cannot be formed for any type of non-linearity.In the context of the cross-diffusion waves for hydro-poro-mechanical coupling (Hu et al., 2020) real stationary solitons, which propagate with a constant velocity without changing their form, are exact solutions of the Korteweg-deVries equation (Regenauer-Lieb et al., 2013a;Veveakis and Regenauer-Lieb, 2015).They depend on the fluid-and solid self-diffusion coefficients only (Veveakis and Regenauer-Lieb, 2015).Cross-diffusion waves are quasi-solitons where the two additional cross-diffusion coefficients cannot be eliminated (Hu et al., 2020).These additional time-dependent properties lead to interesting dynamics.
Following dynamic properties of quasi-solitons have been identified (Zakharov et al., 2004): Quasi-solitons live only for a finite time and can be compared to unstable particles in nuclear physics.Unlike true solitons, quasi-solitons loose energy through their oscillatory tails which can have different wavenumbers in the forward and backward direction of their motion.If the amplitudes of the tails are small, quasi-solitons can be treated as slowly decaying real solitons which lose their energy by radiating quasi-monochromatic waves with wavenumber k 0 in the backward direction.
The discrete particle-like behaviour can be explained by their unusual dispersion relation.Quasi-solitons travel with a constant group velocity ω .When their phase velocity ω(k)/k exhibits a local minimum at a nonzero wavenumber a gap in the spectrum ω(k) appears.According to Zakharov et al. (2004), this peculiar discretization of wave energy has been noticed in many disciplines.Different nomenclatures of quasi-solitons have been adopted.In the theoretical physics community, they have been attributed to Cherenkov radiation.In the ice-wave community, they have been called ice-waves with decaying oscillations and in shallow water theory, they have been identified as capillary-gravity waves.If the amplitudes of the quasi-solitons time/spatial macroscale.This occurs when cross-diffusion waves converge to standing wave solutions when their radiative tails become vanishingly small such as in the non-linear Schrödinger equation or the Korteweg-deVries equation.The approach reveals that instabilities based on the shear and volumetric response of the material at the meso-scale are fundamentally important and have been overlooked.We have shown that incompatibilities of meso-scale accelerations with the overall stress field lead to the nucleation of cross-diffusion waves which travel in an unstable particle-like state with characteristic material velocities c defined by the competition of meso-scale reaction-diffusion processes at the propagating wavefront.The physics of this phenomenon is discussed further in (Regenauer-Lieb et al., 2020).These velocities characterise the progress of internal material timescales for the formation of multiscale space/time dissipative structures and are characteristic properties for the dynamic behaviour of a given material.These internal material clocks are here introduced as a multiscale THMC cascade for coupling the physics of the very small to the very large.
(iii) Cross-diffusion waves have first been discovered for interactions in quantum mechanics such as in photonics where they show anomalous dispersion patterns that, unlike solitons radiate energy in the form of oscillatory (Cherenkov)-tails (Zakharov and Kuznetsov, 1998;Paschotta, 2008).This unusual energy radiation property differentiates them from solitons as quasisolitons.Since they assume an unstable particle-like state (see the example in Fig. 2), the reflections and collisions, when they happen, can lead to a variety of responses (Biktashev and Tsyganov, 2016;Zakharov et al., 2004;Lioubashevski et al., 1996).Despite their nucleation through discrete internal micro-dissipative mechanisms, cross-diffusional waves also show proper soliton wave-like behaviour and can penetrate through each other and reflect from boundaries.However, unlike true solitons, their amplitude and speed are not controlled by initial conditions but by material properties (Tsyganov and Biktashev, 2014).The effect of cross-diffusion is to trigger cross-diffusion waves for critical conditions.They form by THMC-feedback as discrete material instabilities which can be either observed as a local, discrete failure or as damage waves.
They are found to propagate and reflect from boundaries in a multiscale energy cascade which in extreme cases can lead to the formation of a local turbulence phenomenon that seems to appear from nowhere.The phenomenon of wave sampling of energy into ultra-localized events is well known in many disciplines and appears in oceans as 'rogue waves' or in laser physics as sudden 'light-bursts'.In part 2 (Regenauer-Lieb et al., 2020) we will elaborate on the constraints that may lead to an earthquake as a "rogue wave" phenomenon through pumping wave energy over multiple length-scales by multiphysics THMC-feedbacks.
In analogue quantum and hydrodynamic systems the collision, merging and collapse of quasi-solitons (cross-diffusion) waves provides an important mechanism for a direct cascade of energy, which has been also labelled "quasisolitonic turbulence" (Zakharov et al., 2004).This phenomenon may be considered a potential source mechanism for the nucleation of earthquakes.
Their unstable non-local nature and their capability to focus energy from the large-scale to feed local scale 'rogue waves' (Zakharov et al., 2004;Akhmediev et al., 2009) is well documented in photonics (Solli et al., 2007) and shallow water theory (Akhmediev et al., 2009).

Figure 1 .
Figure 1.Dissipative patterns in chemical and biological systems: Simulated and real patterns on Puffer Fish (Sanderson et al., 2006) and strangely patterned siltstone (Zebra Rock, Ranford formation in the Kununurra district in the East Kimberley region of Western Australia).

Figure 2 .
Figure 2. Dissipative patterns in fluid systems: Water molecules exhibit a discrete quantum-like solitary state when forced by a mechanical shaker at a critical condition (here 41 Hz).Periodic finger-like solitary states travel from right to left at a constant velocity.Each snapshotshows 20 ms intervals.Unlike classical solitons their appearance is particle-like.They can pass through each other with a slight loss of amplitude, or 'collide' to create a new state whose direction of propagation is at an angle to that of the original states or disintegrate upon collision (image fromLioubashevski et al. (1996)).

Figure 3 .
Figure 3. Dissipative patterns in solid systems: Two generations of compaction and shear enhanced compaction bands (hand lense for scale) in silt-and sandstones of the Miocene Whakataki formation, south of Castlepoint, North Island, New Zealand.The porosity reduction renders the bands less susceptible to erosion causing a marked positive relief.The periodic bands reminiscent of a standing wave phenomenon are thought to be the result of slow compression of the porous siltstone in the compressive environment (accretionary wedge) of the Hikurangi subduction zone.

Figure 4 .
Figure4.Acceleration waves can originate at a body surface when the existing internal stress gradient is dynamically incompatible with accelerations imposed on particles of the surface.A propagating plane-wave front is shown here for reference, but a plane-wave is not a necessary restriction.Across these surfaces particle accelerations and spatial gradients of velocity are momentarily discontinuous while the velocity itself is continuous.
exchange rate from the A to B phase and vice-versa.In the meso-scale formalism we need to consider information from the local scale processes in the THMC diffusion matrix and decompose the processes leading to the local mass production ξlocal i and ξlocal j .In order to specify this further we define the global volume fraction of the A-phase as: as solitons when the wave amplitude |ψ(x, t)| = |β(x − vt)| propagates without change of form and v and Ω are constants.

Table 1 .
Generalized Thermodynamic Fluxes and Forces in a THMC coupled system (1-D)