A monolithic model for phase-field fracture and waves in solid-fluid media towards earthquakes
AA monolithic model for phase-field fracture and wavesin solid-fluid media towards earthquakes OM ´ A ˇ S R OUB ´ I ˇ CEK
AND R OMAN V ODI ˇ CKA . Abstract
Coupling of rupture processes in solids with waves also propagating in fluids is a prominent phenomenon arisingduring tectonic earthquakes. It is executed here in a single ‘monolithic’ model which can asymptotically capture bothdamageable solids (rocks) and (visco-)elastic fluids (outer core or oceans). Both ruptures on pre-existing lithosphericfaults and a birth of new faults in compact rocks are covered by this model, together with emission and propagationof seismic waves, including, e.g., reflection of S-waves and refraction of P-waves on the solid-fluid interfaces. Arobust, energy conserving, and convergent staggered FEM discretisation is devised. Using a rather simplified variantof such models for rupture, three computational experiments documenting the applicability of this approach arepresented. Some extensions of the model towards more realistic geophysical modelling are outlined, too.
Keywords : Fracture of faults, tectonic earthquake dynamics, elastic waves, elastic-fluid/solid interaction, numer-ical modelling
PACS : 46.50.+a, 91.30.Cd, 91.30.Px
AMS Class. : 74F10, 74J10, 74R20, 74S05, 86-08
Dynamic fracture mechanics is an area of continuum mechanics of solids with wide applications in engineering andparticularly also in geophysics. Global geophysical models typically deal with several very different phenomena andcouple various models due to the layered character of terrestrial planets (including our planet Earth as well as ourMoon). This paper demonstrates the philosophy that a single model can be used instead of several specialized modelsthat otherwise would need to be mutually coupled in a rather complicated way. Such a single universal (we say“monolithic”) model can also be straightforwardly implemented on computers omitting any interfaces which usuallycomplicate implementations. In particular, the transient conditions on fluid/solid interfaces are automatically involvedand need not be specified. This can be a sound advantage of the presented model. Of course, computationally, sucha monolithic-type model may not always make it easier to produce really relevant simulations on the computers wehave at our disposal nowadays. Routine calculations are performed separately either for earthquake sources locally(see, e.g., [10]) or for purely seismic global 3D models, cf., e.g., [25, 26, 29, 46]. Their mutual coupling is alreadytreatable well in the literature [18, 20, 21, 36]. Therefore, there is promising potential that the presented couplingmonolithic approach may become even more amenable in the future with ever increasing computer efficiency.The phenomena we have in mind in this paper involve global seismicity and tectonics . In particular the latterinvolves, e.g., ruptures of lithospheric faults or a birth of new faults, generating seismic waves which then propagatethrough the solid-like silicate mantle and iron-nickel inner core both in the shear (S) or the pressure (P) modes. Incontrast to the P-waves (also called primary or compressional waves), the S-waves (also called secondary waves) aresuppressed in the fluidic iron-nickel outer core and also in the water oceans.A very low attenuation of seismic waves is the ultimate phenomenon. We, therefore, take the
Jeffrey’s rheology model (i.e., a serial combination of the Maxwell and Kelvin-Voigt rheology as in Figure 1/bottom-left, cf., e.g., [33]or also [39]) as a basic global “monolithic” ansatz. In various limits in the deviatoric and the volumetric parts, wemodel different parts of planet Earth. Jeffrey’s rheology also seems more realistic in particular because it covers (inthe limit) also the Kelvin-Voigt model applied to the volumetric strain (actually considered as a starting model inFigure 1/left) whereas the pure Maxwell rheology, allowing for big creep during long geological periods, is not arelevant effect in the volumetric part. The support from the grants 17-04301S (as far as dissipative evolutionary aspects concerns) and 19-04956S (as far as modelling of dynamicand nonlinear behaviour concerns) of the Czech Sci. Foundation and VEGA 1/0078/16 of the Ministry of Education, Science, Research and Sportof the Slovak Republic, and the institutional support RVO:61388998 ( ˇCR) are acknowledged. Institute of Thermomechanics, Czech Academy of Sciences, Dolejˇskova 5, CZ-182 00 Praha 8, Czech Republic. Mathematical Institute, Charles University, Sokolovsk´a 83, CZ-186 75 Praha 8, Czech Republic Civil Engineering Faculty, Technical University of Koˇsice, Vysokoˇskolsk´a 4, SK-042 00 Koˇsice, Slovakia a r X i v : . [ phy s i c s . g e o - ph ] A ug especting the solid parts of the model, we use the Lagrangian description even in the fluid regions, i.e. hereall equations are formulated in terms of displacements rather than velocities. The reference and the actual spaceconfigurations automatically coincide with each other in our small strain (and small displacement) ansatz, which iswell relevant in geophysical short-time scales of seismic events.In the solid-like part, various inelastic processes are considered to model tectonic earthquakes on lithosphericfaults together with long-lasting healing periods in between them, as well as aseismic slips, and various other phe-nomena. To this goal, many internal variables may be involved such as aging/damage, inelastic strain, porosity, watercontent, breakage (i.e., essentially like another damage-like internal variable in modeling of granular materials), andtemperature, cf. [30, 31, 33]. On the other hand, those sophisticated models are focused rather on local events aroundthe tectonic faults without ambitions to be directly coupled with the global seismicity. Here, rather for the lucidityof the exposition, we reduce the set of internal variables to only one scalar variable and one matrix-valued variable,namely damage/aging and an inelastic strain, respectively. The calculations in Sect. 5 involve the only former one.Even this simple scenario, however, has a capacity to trigger a spontaneous rupture (so-called dynamic triggering )with emission of seismic waves and, in a certain simplification, can serve as a seismic source coupled with the overallglobal model. The mentioned inelastic strain can capture Maxwell-type rheologies relevant in the solid mantle andinner core to capture long-term creep (aseismic) effects up to 10 yrs.Let us emphasize that the usual models are focused only on either the propagation of seismic waves along thewhole globe while their source is considered given, or on the description of seismic sources due to tectonic events,but not their mutual coupling. If a coupling is considered, then it concerns rather local models not considering thelayered structure of the whole planet, cf. e.g. [5,6,20,21,32]. The reality ultimately captures very different mechanicalproperties of different layers of the Earth, in particular the mantle and the inner core which are solid on the short-timescales versus the outer core and the oceans which are fluidic even on the short-time scales. Some other coupledmodels implement a pre-existing fault (in contrast to our model allowing for new faults birth, cf. Section 5.2) and usea slip-dependent friction law [45].The goal of this contribution is threefold: α ) to present a model that might simultaneously capture the sources of seismic waves that necessarily behavenonlinearly (like ruptures of tectonic faults) and the propagation of seismic waves possibly even over the wholeplanet (when further special algorighmic and computational methods would be used, not handled in this article,however), both phenomena being mutually coupled. β ) by proper scaling to approximate viscoelastic (so-called Boger’s [7]) or merely elastic fluids that are relevant inthe outer core and in the oceans (with a very low or just zero viscosity) where S-waves cannot propagate while P-waves are only refracted on the solid/fluid interfaces (in particular on Guttenberg’s core-mantle discontinuity). γ ) to document computational efficiency of the monolithic model at least on 2-dimensional rather local simulationsof various, quite distinct geophysical events.We refer to [40] for the rigorous analysis as far as asymptotics of the scaling of viscosities in even a more complexmodel involving also self-induced gravity, tidal, Coriolis, and centrifugal forces; cf. Sect. 6.2. Mechanical models in general (and those used in geophysics in particular) typically are (or should be) believed to begoverned by energies and, most often, in a way that the conservative and the dissipative parts are separated. In theisothermal variant, the systems have a simple general structure of an abstract dissipative dynamic equation M (cid:48) .. q + R (cid:48) . q ( q , . q ) + E (cid:48) ( q ) = F ( t ) (1)with a kinetic energy M , a (pseudo) potential of dissipative forces R ( q , · ) , a stored energy E , and external forcing F as a time-dependent functions of the state q in the reference domain Ω . In (1), R (cid:48) . q denotes the partial differentialof R = R ( q , . q ) with respect to . q . In fact, R ( q , · ) may be non-differentiable at . q = R (cid:48) . q ( q , · ) may be set-valued and (1) is to be an inclusion rather than equality;yet we will ignore these technicalities in this presentation.2his state q typically involves, beside of the displacement, also some internal variables like the inelastic strainsdescribing creep and a “permanent” deformation resulting from re-occurring shifts during earthquakes, damage/aging,etc., complying with the concept of generalized standard materials with internal variables [15]. In (1), we use thenotation . q = dd t q and .. q = d d t q , and ( · ) (cid:48) denotes the differential. Also, (1) takes the structure of dissipative Hamiltoniansystem , and the Hamilton variational principle [17] extended to dissipative systems as in [3] says that the solution q to (1) on a fixed time interval [ , T ] is a critical point of the integral functional q (cid:55)→ (cid:90) T M ( . q ) − E ( q ) − (cid:10) f , q (cid:11) d t (2)with a nonconservative force f = R (cid:48) . q ( q , . q ) − F ( t ) considered fixed on the affine manifold respecting some initial orterminal conditions. The notation (cid:104)· , ·(cid:105) in (2) means the value of a functional (in our case f ) on a test function (in ourcase q ), which for sufficiently smooth functions can be understood as an integral, see e.g. (9d) below.This energy-governed structure (1) allows to control the energetics: indeed, testing (1) by . q and using the calculus (cid:104) M (cid:48) .. q , . q (cid:105) = dd t M ( . q ) and (cid:104) E (cid:48) ( q ) , . q (cid:105) = dd t E ( q ) , we arrive (at least formally) to the energy balance on a time interval [ , t ] : M ( . q ( t ))+ E ( q ( t )) (cid:124) (cid:123)(cid:122) (cid:125) kinetic + stored energyat time t + (cid:90) t Ξ ( q ( s ) , . q ( s )) d s (cid:124) (cid:123)(cid:122) (cid:125) dissipated energy overthe time interval [ , t ] = M ( . q ( ))+ E ( q ( )) (cid:124) (cid:123)(cid:122) (cid:125) kinetic+stored energyat time 0 + (cid:90) t (cid:10) F ( s ) , . q ( s ) (cid:11) d s (cid:124) (cid:123)(cid:122) (cid:125) work done by loadingover time interval [ , t ] (3a)with the dissipation rate Ξ (cid:0) q , . q (cid:1) = (cid:10) R (cid:48) . q ( q , . q ) , . q (cid:11) . (3b)Moreover, this structure allows for various numerically stable time discretisations, and for rigorous analysis as faras convergence and existence of solutions to (1) concerns. Denoting the displacement by u , we use the usual concept of small strains, with the total strain e = e ( u ) = ∇ u + ( ∇ u ) (cid:62) and the inelastic, symmetric, trace-free strain π . As standard in small-strain theories, we consider the additivesplit, sometimes called Green-Naghdi split [13], of the total strain e = e ( u ) as e ( u ) = e el + π (4)with e el denoting the elastic strain and π the inelastic strain, the latter being considered as trace free. To distinguishthe response on the shear or compression load, we further introduce its orthogonal decomposition to the spherical andthe deviatoric part: e ( u ) = sph e ( u ) + dev e ( u ) with dev e = e − tr e I , (5a)dev e = dev ( e el + π ) = e dev + π with e dev = dev e el , (5b)where I is the identity matrix while tr e = ∑ i = e ii denotes the trace of e . A rather minimal scenario for a simplified seismic source model in the solid part (the mantle) is an isothermal damage(also called aging) in the deviatoric part. To model also the aseismic slip (creep) or a permanent displacement duringcumulation of subsequent shift (which is important in particular if healing is considered in between re-occurringearthquakes), as well as a low attenuation of seismic waves, the Maxwell viscoelastic rheology is considered in thedeviatoric response. For well-posedness of the problem, it is suitable to combine it with the Kelvin-Voigt rheology,which results to the Jeffreys rheological model, as used in [33]. The spherical response in the solid part (the mantle3 ρ PSfrag replacements e ( u ) e ( u ) π αασ σ sph σ sph σ dev σ dev σ dev ≡ e ( u ) dev e ( u ) e dev dev e el sph e ( u ) sph e ( u )sph e el G E ( α ) G MX G V G V ( α ) G M K E K E K MX K V K V DAMAGEABLE VISCOELASTIC SOLID VISCO ELASTIC(BOGER’S) FLUIDSOLID / FLUIDINTERFACE DEVIATORICPART OFTHE MODELSPHERICAL(VOLUMETRIC)PART OFTHE MODELUNIVERSALLY USABLE PART
Figure 1:
Schematic diagram for solid-fluid interaction. The spherical response is the Kelvin-Voigt rheology both inthe solid and in the fluid parts. The deviatoric response of the solid uses a viscoelastic Jeffreys rheology subjected todamage α , while the fluid uses the mere Stokes model. and inner core) exhibits just the undamageable Kelvin-Voigt rheology, reflecting the phenomenon that neither any bigpermanent deformation nor damage by compression is possible. Cf. Figure 1/left for the overall rheological model.The state q = ( u , π , α ) then involves the displacement u , the (trace-free) inelastic strain π , and the scalar damagevariable α ; alternatively, α is called also aging and we follow the usual convention that α = α = specific stored energy is postulated as ϕ ( e , π , α ) = K E | sph e | + G E ( α ) | e dev | + γ DAM ( α ) = K E | sph e | + G E ( α ) | dev e − π | + γ DAM ( α ) (6)with K E the elastic bulk modulus and G E = G E ( α ) the elastic shear modulus introduced for a damageable Lam´e mate-rial, and γ DAM = γ DAM ( α ) damage stored energy (also understood as the damage toughness). The specific dissipationpotential is ζ ( α ; . e , . π , . α ) = ζ VSC ( α ; . e , . π ) + ζ DAM ( . α ) with (7) ζ VSC ( α ; . e , . π ) = K V | sph . e | + G V ( α ) | dev . e − . π | + G M | . π | and with ζ DAM = ζ DAM ( . α ) ≥ ζ VSC .The former one is the potential of non-conservative forces related to damage evolution while the latter one togetherwith elastic parts in the stored energy introduces Kelvin-Voigt rheology for the spherical part and Jeffrey rheologyfor the deviatoric part. We consider still the specific kinetic energy in the usual form ρ | . u | with ρ = ρ ( x ) the massdensity.To introduce a length-scale into the model that allows for “controlling” a typical width of the damaged zonearound fault core (usually varying in between 1–100 m) and of the narrow cataclasite fault core (i.e. the principalslip zone, usually varying in between 0.1–10 m) where the inelastic strain is accommodated, the gradient theories areused. This means that ϕ from (6) is augmented both by | ∇α | and | ∇π | terms. In a special choice γ DAM ( α ) = G c 12 ε α together with G E ( α ) = ( ε / ε +( − α ) ) G where G is the shear modulus of the undamaged material (up to a smallvalue of ε / ε ) and ε > ε ), the “augmented” specific stored energy ϕ A takes the form ϕ A ( e , π , α , ∇π , ∇α ) : = K E | sph e ( u ) | + (cid:16) ε ε +( − α ) (cid:17) G | dev e ( u ) − π | + κ | ∇π | + G c (cid:16) ε α + ε | ∇α | (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) crack surface density (8)with ε > κ (dimension [J/m]) allows to control the width of the mentioned cataclasite fault corewhile ε (dimension [ m ] ) allows to control the width of the mentioned damaged zone around the fault core, and aso-called fracture toughness G c > ]) fixed. This is, in fact, the usual phase-field crack model,4ere considered only in the deviatoric component reflecting the phenomenon that pure tension (and Mode I cracks)is avoided in geophysically relevant situations while compression does not lead directly to cracks, so that only shearmay lead to rupture (in Mode II). In the static situation, for ε converging to zero, this so-called Ambrosio-Tortorellifunctional approximates the Griffith crack model [14], originally shown in the scalar case [1] and later for the vectorialcase by [12]. In the dynamical case, this approximation property is not rigorously justified, but nevertheless thisphase-field model is routinely used, cf. e.g. [8, 28, 35, 41]. Such a crack model is relevant particularly if a new fault isnucleated in the compact rock, although it is a rare event geophysically.
Besides solid regions, there are layers in our planet that are naturally fluidic, namely the outer core and, of course, theoceans. From the seismic wave propagation, they must exhibit some spherical elastic response (otherwise, ideally,incompressible fluids would lead to unphysically infinitely large P-wave speed) but no deviatoric elastic response toprevent S-wave propagation in such regions. In any case, damage becomes irrelevant in these regions.One scenario is to keep still a (presumably small) viscous response both in the spherical and the deviatoric part,leading respectively to a combination of the Kelvin-Voigt solid rheology and the Stokes fluid rheology. In the notationfrom Section 3.1, this means that we take G E = γ DAM = G V > α , κ = G c =
0, and we formallyput G M = ∞ . The last action results to π = π | t = =
0. In fact, an equivalent optionwould be to put G V = ∞ and keep G M > viscoelastic ( Boger ’s [7]) fluid ,cf. Figure 1/right-part. The S-waves actually can slightly penetrate into such fluids but are soon attenuated and thuscannot propagate practically.Another scenario is to suppress any viscosity, obtaining, thus, the merely elastic fluid . This can be achieved, inaddition to taking G E = γ DAM = κ =
0, and G c =
0, by putting both K V = G V =
0. Such fluids are fullyconservative, leading to a linear hyperbolic problem. The S-waves cannot propagate through such fluidic regions atall.
Let us consider the reference domain Ω occupied by a viscoelastic body in question (e.g., the whole planet Earth).Of course, the local potentials ϕ A and ζ , and also the mass density ρ are allowed (and supposed) to depend on x ∈ Ω ,and in particular, they can vary on the solid and the fluid domains; let us denote them by Ω S and Ω F , respectively. Wewill not indicate this dependence explicitly for notational simplicity.Taking into account the additive split (4), the overall functionals used in Section 2 are now: E ( q ) = E ( u , π , α ) = (cid:90) Ω ϕ A ( e ( u ) , π , α , ∇π , ∇α ) d x , (9a) R ( q , . q ) = R ( α ; . u , . π , . α ) = (cid:90) Ω ζ VSC ( α ; e ( . u ) , . π ) + ζ DAM ( . α ) d x , (9b) M ( . q ) = M ( . u ) = (cid:90) Ω ρ | . u | d x , (9c) (cid:104) F ( t ) , u (cid:105) = (cid:90) Ω g ( t ) · u d x (9d)where g is a bulk force, typically the gravitational force independent of time or some more general force as, e.g.,in (23) or (24) below. Note the (9d) does not involve the internal variables ( π , α ) , in accord with the conventionalconcept that the internal variables are not directly subjected to outer forcing.More specifically, the system that results from (1) by considering (9) together with (8) looks as: ρ .. u − div ( σ sph + σ dev ) = g in Ω (10a)with σ sph = K V sph e ( . u ) + K E sph e ( u ) in Ω , (10b)and σ dev = G V ( α ) (cid:0) dev e ( . u ) − . π (cid:1) + G E ( α )( dev e ( u ) − π ) in Ω S , G V dev e ( . u ) in Ω F , (10c) G M . π = σ dev + div ( κ∇π ) in Ω S , (10d)5 (cid:48) DAM ( . α ) + G (cid:48) E ( α ) | ( dev e ( u ) − π ) | + G c ε α = div ( ε G c ∇α ) in Ω S . (10e)Note that the spherical and the deviatoric stresses are orthogonal, i.e. σ sph : σ dev =
0. Also note that π and α arerelevant only on the solid domain Ω S , cf. also Figure 1.There are no transient conditions for displacement/stress on the solid-fluid interior interfaces (i.e. between outercore and mantle and inner core and possibly also between mantle and oceans) because (10a) is considered on the wholedomain Ω . In particular, the stress-vector equilibrium and continuity of the displacement across these interfaces areautomatically involved and does not need to be written explicitly. On the other hand, there are boundary conditionsto be read from the abstract equation (1), namely ( σ sph + σ dev ) (cid:126) n = ∂Ω , (11a) κ ∂π∂ (cid:126) n = ε G c ∂α∂ (cid:126) n = ∂Ω S , (11b)where (cid:126) n denotes the normal to the boundary ∂Ω of Ω (e.g. the surface of the Earth) the interior boundaries ∂Ω S (i.e.on the mantle/core and mantle/oceans and inner/outer-core interfaces) and ∂∂ (cid:126) n = (cid:126) n · ∇ denotes normal derivative.We have in mind an initial-value problem, so that we have to complete the system still by initial conditions: u | t = = u and . u | t = = v in Ω , (12a) π | t = = π and α | t = = α in Ω S . (12b)Although it is usually not an aspect under attention in geophysical modelling, let us mention that it is possible toprove rigorously that, under a suitable data qualification, the initial-boundary-value problem (10)–(12) has a solutionwhich also satisfies the energy conservation (3), cf. [40] for the case G V independent of α while [27, Sect. 7.5] outlinesmodifications if G V ( · ) = τ V G E ( · ) for some relaxation time τ V = τ V ( x ) > One can approximate the fluid models from Sect. 3.2 asymptotically when sending corresponding parameters of thesolid model in Sect. 3.1 to their limits. More specifically, Boger’s viscoelastic fluid can be approached from theJeffreys solid by sending G M → ∞ , G E → , κ → , and ζ DAM → Ω F . (13)When π = π converges to 0 on the fluidic domain Ω F . Under suitable data qualificationand scaling, this convergence can rigorously be justified together with convergence of the energy balance, cf. [40].Sending further G V → K V → Ω F , (14)we approach the merely elastic fluid in Ω F . Again, this convergence can rigorously be proved together with conver-gence of the energy balance (i.e. the energy dissipated via viscous attenuation in the Bogger fluid actually convergesto zero) but for a slightly modified model with the bounded elastic stress and constant G V in the solid part Ω S , cf. [40].The idea behind the “monolithic” approach is to take, instead of the limit fluid in Sect. 3.2, its approximation andto implement computationally only the solid model in the whole domain Ω . Again, the interface conditions between Ω S and Ω F expressing here continuity of displacements and the traction stresses will thus be covered automaticallywithout any extra effort in coding. We write the abstract equation (1) as a first-order system, which is more suitable for time discretisation than theoriginal 2nd-order system, in particular because it allows for varying time steps during simulations. In view of (9), ithas the structure . u = v , (15a)6 (cid:48) . v + R (cid:48) v ( α ; v , . π ) + E (cid:48) u ( u , π , α ) = F ( t ) , (15b) R (cid:48) . π ( α ; v , . π ) + E (cid:48) π ( u , π , α ) = , (15c) R (cid:48) . α ( . α ) + E (cid:48) α ( u , π , α ) = . (15d)It is important that R is additively split, cf. (9b), which suggests the splitting with the time discretisation of the state q to the components ( u , π ) and α and in the same manner also the calculation is performed. Although E is necessarilynonconvex to facilitate modelling of sudden rupture events, it is also advantageous that both E ( · , · , α ) and E ( u , π , · ) are convex or, here in the ansatz (8), even quadratic.Considering a time step τ = ( τ k ) k ∈ N with k = , , .... indexing time levels t k in the time-discrete system so that τ k = t k − t k − , and using also the Crank-Nicholson mid-point strategy, we devise the staggered (also called fractional-step splitting) system u k τ − u k − τ τ k = v k − / τ with v k − / τ = v k τ + v k − τ , (16a) M (cid:48) v k τ − v k − τ τ k + R (cid:48) v (cid:16) α k − τ ; v k − / τ , π k τ − π k − τ τ k (cid:17) + E (cid:48) u ( u k − / τ , π k − / τ , α k − τ ) = F ( t k ) , (16b) R (cid:48) . π (cid:16) α k − τ ; v k − / τ , π k τ − π k − τ τ k (cid:17) + E (cid:48) π ( u k − / τ , π k − / τ , α k − τ ) = , (16c) R (cid:48) . α (cid:16) α k τ − α k − τ τ k (cid:17) + E (cid:48) α ( u k τ , π k τ , α k − / τ ) = . (16d)The system (16) is to be solved, recursively, for k = , , ... , starting from the initial conditions u τ = u , v τ = v , π τ = π , and α τ = α , cf. (12). The system is decoupled in the sense that the calculation is performed separately for ( u k τ , v k τ , π k τ ) from (16a-c) and α k τ from (16d), as mentioned already above.This scheme exhibits energy conservation [42], which can be seen by multiplying (16b) by the mid-point velocity v k − / τ , (16c) by the inelastic-strain rate ( π k τ − π k − τ ) / τ k , and (16d) by ( α k τ − α k − τ ) / τ k . For the M -term in the formertest, we use (16a) and the binomial formula to obtain an analog of the calculus (cid:104) M (cid:48) . v , v (cid:105) = dd t M ( v ) as the equality (cid:68) M (cid:48) v k τ − v k − τ τ k , v k − / τ (cid:69) = (cid:68) M (cid:48) v k τ − v k − τ τ k , v k τ + v k − τ (cid:69) = M (cid:0) v k τ (cid:1) − M (cid:0) v k − τ (cid:1) τ k , (17)while another binomial formula for the quadratic functional E ( · , · , α k − τ ) gives (cid:68) E (cid:48) ( u , π ) ( u k − / τ , π k − / τ , α k − τ ) , (cid:16) v k − / τ , π k τ − π k − τ τ k (cid:17)(cid:69) = E ( u k τ , π k τ , α k − τ ) − E ( u k − τ , π k − τ , α k − τ ) . (18)Similarly, the latter test gives (cid:68) E (cid:48) α ( u k τ , π k τ , α k − / τ ) , α k τ − α k − τ τ k (cid:69) = E ( u k τ , π k τ , α k τ ) − E ( u k τ , π k τ , α k − τ ) . (19)Summing (18) with (19), we can enjoy cancellation of ± E ( u k τ , π k τ , α k − τ ) . Then, summing it over k = , , ... , weobtain a discrete analog of the energy balance (3) as an (exact!) equality, not only as an inequality. This eliminatesa spurious numerical attenuation usually exhibiting by implicit discretization schemes and is helpful during codingbecause accidental mistakes can thus be immediately detected.After applying still a space discretisation, the recursive scheme (16) can be implemented on computers, leadingto one linear and one linear-quadratic-programming problem at each time level. Therefore, this scheme can be solvedby finite (non-iterative) algorithms, is robust (i.e. numerically stable), convergent, energy conserving. It can beinterpreted as a combination of the (generalized) Newmark time discretisation with staggered (fractional-step split)for the damage flow-rule, cf. also [19] for a similar scheme. Here, for the space disretisation in the examples inthe following Section 5, we used just the simplest P1 finite elements. On the other hand, a more detailed numericalanalysis as far as rate of convergence or error estimates is another challenge not addressed in this article, and it isstandardly considered as very difficult in such highly nonlinear problems, however.7 Sfrag replacements k m
200 km
PSfrag replacements
100 km200 km
Figure 2:
An example of triangular meshing of the rectangular domain used for Figure 3 below. The element size isabout . We present rather academic (not entirely in real geophysical scale and only 2-dimensional) computational examplesdocumenting efficiency and applicability range of the above presented monolithic model even in its simplified variantwhen G M → ∞ so that π →
0. In other words, we neglect the inelastic strain, which is relevant for short-time/rangeevents like ongoing ruptures and earthquakes around hypocentres if the (low) Maxwellian attenuation or aseismiccreep are neglected. Also, for the illustrative calculations, we have neglected the viscosity, which is justified in thefluidic part (where it leads to the linear hyperbolic problem for a merely elastic fluid [40]) but not in the nonlinearsolid part. Regardless, there is a belief that this shortcut does not influence substantially the presented simulations.We consider a solid rock with the elastic bulk modulus K ES =
600 GPa, the shear modulus G ES =
250 GPa, anddensity ρ S = − . Thus, in this two-dimensional situation, the speed of the P-waves (= the sound speed) is v PS = (cid:112) ( K ES + G ES ) / ρ S = .
04 kms − , and the speed of the S-waves is v SS = (cid:112) G ES / ρ S = .
07 kms − . As for thefluids, we consider the elastic bulk modulus K EF = ρ F = − which provides thespeed of the P-waves v PF = (cid:112) K EF / ρ F =
10 kms − . The other parameters of the model which influence the damagepropagation are G c = -2 , ε = − m, ε = v g = − .Moreover, the coefficient 3 in (6)–(8) and (10b) is to be 2 in the two-dimensional setting. Also in these three calcula-tions, the interface between solid and fluid is considered distinctly, with the intent to provide various possibilities forreflected/refracted waves. Nevertheless, it should be emphasized that in all three calculations the same computationalmodel is used and the differences appear only in initial/boundary conditions and geometry of the solid and the fluidphases.The space discretisation of the 2-dimensional body made by the P1 finite elements includes an irregular but moreor less uniformly sized mesh, whose smallest element size h min = h min is required at the interface between solid and fluid and also on pre-existingfaults. The largest elements close to the top face are only a half larger. The size of the mesh (about 2 km) was chosenin accordance with computational hardware (namely Intel(R) Core(TM) i5 CPU 2.3GHz and 8GB RAM) used tocompromise the visualization of the expected geophysical events with CPU time and memory needed. Although, theshown mesh pertains to the first calculations. the second and the third ones are very similar with the same smallestelement size h min so that they are not shown explicitly.The time discretisation (16) is implemented with adaptively varying time step τ k shrinking when fast rupture andwave propagation start. To capture the latter phenomenon properly, the so-called Courant-Friedrich-Lewy (CFL)condition [11] has been respected, i.e. here τ k < h min / max ( v PS , v PF ) ; in fact, the CFL-condition is ultimately needed8 (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) PSfrag replacements k m
200 km Ω F Ω S partly pre-damagednarrow stripe - a faultepicentrehypocentresolid-fluidinterface D A M A G E A B L E S O L I D E L A S T I C F L U I D s o li d / fl u i d i n t e r f a c e DEVIATORICPART OFTHE MODELSPHERICAL(VOLUMETRIC)PART OFTHE MODELUNIVERSALLYUSABLE PART
Figure 3:
A computational 2-dimensional domain and boundary conditions for the experiment in Sect. 5.1. Uppertrapezoidal part is the solid from Figure 1/left-part while the lower part is the fluid from Figure 1/right-part. The uppersolid part is shifted by the boundary conditions to the left while the lower solid/fluid part is horizontally constrained.At time t = , the fault is “compact” only at the middle part (where the hypocentre of an earthquake will be) whilethe rest of the fault (i.e. both sides) is (partly) damaged. for explicit time discretisation to ensure stability and convergence but also desired for our implicit one when wavesare emitted and propagating. In particular, we took τ k =
100 s for the initial phase of the prescribed loading and thendecreased to τ k = .
146 s when the rupture and wave propagation were triggered. In absolute values, the time instant t of the switch to the shorter time step was chosen as follows: the first calculation t =
32 ks, the second calculation t =
49 ks, the third calculation t = . × in the performed academic calculations. In real structures it can be severalorders higher. The prescribed horizontal shift of the upper plate at the moment of rupture initiation regarding thevelocity v g is in any case more then 27 m which is almost 10 smaller then the size of the considered domain. The most typical scenario leading to tectonic earthquakes is that a fault (as a flat usually straight stripe with weakeneddamage threshold) is gradually stretched until the shear stress reaches the critical value to trigger the rupture. Thena fast shift occurs, within which a seismic (mainly S-) wave can be emitted. The geometry of the 2-dimensionalcomputation region with the (a bit hypothetically) horizontally positioned fault together with boundary conditions isdepicted in Figure 3; the boundary conditions imposes a prescribed displacement applied from both sides of the topsolid layer.The results of this example are depicted in Figure 4 in eight selected snapshots after rupture triggering, showingthe spatial distribution of the kinetic and the stored energy in its decomposition to the shear and the spherical parts.The instants are expressed only by increments with respect to the first one t due to the large time scale of the wholecalculation. The actual number of numerical time steps can be computed by doing the increment by τ k = . t and t . The interesting moment is when the upper plate is movedsufficiently far towards left and thus the fault is stretched enough so that the rupture occurs in the middle; i.e. theearthquake starts in the hypocentre. In fact, the fracture toughness G c is put higher close to the boundary Γ to preventrupture starting from the boundary where there is necessarily a stress concentration due to the fast varying boundaryconditions.At that fault rupture moment, the strain energy around is relaxed and a seismic wave is emitted and starts propa-gating through the solid part. One can see that this wave is not a ball-shaped partly because the seismic source (therupturing area) is rather a surface than a point and partly (or mainly) because the seismic source by a slip of the faultgenerates rather an S-wave towards normal (i.e. here vertical) directions (which are slower), while the horizontal-like fronts are rather P-waves (propagating faster). This is also shown in Fig. 5-left. The pertinent waves in kineticenergies around the hypocentre are shown for the snapshot corresponding to the time t . The P-wave is sufficientlyseparated from the S-wave. This is also shown in Fig. 5-left. The pertinent waves in kinetic energies around thehypocentre are shown for the snapshot corresponding to the time t , the detail reveals the S-wave which can also beseen for the same time instant in the shear energy in Fig. 4-middle. The P-wave is well separated from the S-wave.Here we advantageously exploit that in computer we have at disposal the strain energy split into the spherical and9 INETIC ENERGY SHEAR STORED ENERGY SPHERICAL STORED ENERGY t = . s t = t + . s t = t + . s t = t + . s t = t + . s t = t + . s t = t + . s t = t + . s Figure 4:
Simulations of a rupture of a fault in the middle part, emission of a seismic (mainly S-) wave in thehypocentre, its propagation and reflection on the solid/fluid interface where a weak P-wave can be observed in thefluid part, and eventually S-wave reaches also the surface (in an so-called epicentre) where the earthquake manifests.The displacement is magnified 100 × to visualize the deformation. t Figure 5:
Various types of seismic waves recognizable in the kinetic energy drawings for the relevant parts from twoselected snapshots from Fig. 4-left. The specification is possible due to the shear energy (Fig. 4-middle) where onlyS-waves can be visible, while not the P-waves. (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)
PSfrag replacements k m
200 kmΩ F Ω S partly pre-damagednarrow stripe - a faultsolid-fluidinterface D A M A G E A B L E S O L I D E L A S T I C F L U I D solid / fluidinterface Figure 6:
A computational 2-dimensional domain and boundary conditions for the experiment in Sect. 5.2. At time t = , the fault is “compact” at the left part while the right part is substantially damaged. the shear parts, from which one can distinguish these types of waves; in particular, the shear energy in Fig. 4-middleshows clearly the S-wave.When the seismic wave reaches the Earth surface, it is reflected and partly starts propagating to the sides, i.e. theearthquake starts at the epicentre. On the other side, when reaching the solid/fluid interface below the hypocentre,the wave is mainly reflected (which documents that it is rather an S-wave) and to a small extent it generates a slightP-wave in the fluidic part. Some other reflection is seen on the fault itself, which is actually a well recognizedphenomenon related to (or serving for identifying of) the width of the low-velocity damaged fault zone, cf. e.g. [4].The detail plot in Fig. 5-right specifies the form of the reflected waves for the snapshot corresponding to the time t comparing it at least with the shear energy in Fig. 4-middle at the same time instant. Another interesting event, although very rare from the mankind time scale, is a nucleation of damage in compactrocks, giving rise to a new damage regions, i.e. to a new fault. The geometry of the 2-dimensional computation regiontogether with the boundary conditions imposing (unlike the previous example) the increasing displacement only atthe right face of the top solid layer is depicted in Figure 6. It differs from the previous example also by the initialconditions which makes the fault largely damaged on its right-hand part. Therefore, the upper plate can quite easilyslide onto the lower plate while being rather well connected with it on the left-hand part.The results of this example are depicted in Figure 7 in seven selected snapshots of spatial distribution similar to theprevious example. When the upper plate is stretched enough, it starts rupturing on an a-priori not pre-defined place.The energetics of the model dictates that the new damaging area is a (relatively) narrow plane which is positioned atabout 45 degrees with respect to the existing fault. Such position, together with the slip orientation, is referred to asa normal fault , in contrast to reverse (thrust) faults, or strike or vertical faults. In fact, the new fault plane slightlycurves, the dip being steeper near the surface while shallower with increased depth. This is referred to as a listric ormal fault .The rupture starts on the existing fault due to a slight stress concentration and then propagates towards the Earth’ssurface with an increasing speed, emitting a seismic wave. This wave propagates mainly through the upper plate (asan S-wave which can be observed in the plots of kinetic and shear stored energy for the time instants t and t ) butpartly penetrated through the existing horizontal fault into the lower plate and then even towards the fluidic layerwhere it creates a relatively strong P-wave, cf. Figure 7-left the time instant t . Another interesting situation, which may occur together with a dip-slip fault as shown in the previous example, isarising of a reverse (thrust) fault. The geometry of the 2-dimensional computation region together with the boundaryconditions imposing the increasing displacement at both lateral faces of the top solid layer but in opposite directionis depicted in Figure 8. It differs from the previous example also by the initial damage conditions where both endsof the fault are largely damaged. Therefore, the upper plate easily slides on the lower plate at its side parts while thecentral part is well connected with the lower plate and allow for stress concentration.The results of this example are shown in Figure 9 in eight selected snapshots of spatial distribution similarly asin the previous examples. When the upper plate is substantially compressed, it starts rupturing on an a-priori not pre-defined place. The energetics of the model dictates that the new damaging area is a (relatively) narrow plane whichis positioned in about 60 degrees with respect to the existing fault. Such position, together with the slip orientation,is referred to as a reverse (thrust) fault , here again a bit curved (listric) like in the previous normal-fault case. Duringthis first rupture, the newly born thrust fault starts on the existing horizontal fault and then propagates towards theEarth surface with an increasing speed, emitting a seismic wave. This wave propagates mainly through the upperplate but partly penetrates through the existing horizontal fault into the lower plate and then even towards the fluidiclayer where it creates a relatively strong P-wave, cf. Figure 9-left.After a certain time (determined rather by evolving boundary conditions), it causes the horizontal fault rupturewith emission of another S-waves which propagates/reflects in a similar way as in the first example in Sect. 5.1.
The mentioned academical character of the presented, rather simple model can be suppressed when enhancing it byvarious ways.
It is a well recognized effect that the Griffith fracture model, which is approximated by the ansatz (8) if ε →
0, isrealistic as far as fracture (cracks) propagation but has difficulties with nucleation of cracks.This effect can be seen when calculating the damage driving force ϕ (cid:48) α used in (8), obtaining G (cid:48) E ( α ) | e dev | + γ (cid:48) DAM ( α ) = ( − α ) G E0 | e dev | + G c α / ε . Counting with the shear stress σ dev : = dev ϕ (cid:48) e = G E ( α ) e dev , one can see thatthe stress needed for triggering damage (rupture) indeed growths like 1 / √ ε for ε →
0, which causes the mentioneddrawback of the limit Griffith-type model.One way to get rid of this drawback is to consider, instead of G E ( α ) = ( ε / ε +( − α ) ) G E1 used in (8), a moregeneral convex decreasing nonlinearity G E : [ , ] → ( , ∞ ) depending on ε used in the crack surface density in (8),cf. e.g. [9]. When requiring G (cid:48) E ( ) to blow up like 1 / ε , one can still consider ε > G (cid:48) E ( ) and G c . Such softening may lead to localization of damage even without considering the Ambrosio-Tortorellimodel, as shown by [2].Another worthy modification which facilitates nucleation of rupture is to make e el = e − π (cid:55)→ ϕ ( e , π , α ) nonconvexif α is close to 1, reflecting the lost of stability of the rock when damage is sufficiently developed. The 2-homogeneousansatz suggested by [34] and used often in geophysical models is ϕ ( e , π , α ) = λ E | sph e | + G E ( α ) | e − π | − αγ R | e − π | sph e (20)with λ E > G E ( α ) = G − α G R with undamaged shear modulus G like in (8) and 0 < G R < G sensitivity of the shear modulus to damage,12 INETIC ENERGY DAMAGE SHEAR STORED ENERGY t = . s t = t + . s t = t + . s t = t + . s t = t + . s t = t + . s t = t + . s Figure 7:
Simulations of a birth of a new fault in a normal position and a depression on the Earth surface, togetherwith emission of a seismic (mainly S-) wave during completion of the rupture, its propagation and creation of a ratherstrong P-wave in the fluidic domain. The displacement magnified 100 × . (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) PSfrag replacements k m
200 kmΩ F Ω S partly pre-damagednarrow stripe - a faultpartly pre-damagednarrow stripe - a faultsolid-fluidinterface D A M A G E A B L E S O L I D E L A S T I C F L U I D solid / fluidinterface Figure 8:
A computational 2-dimensional domain and boundary conditions for the experiment in Sect. 5.3. At time t = , the fault is “compact” at the central part only while both the left and the right parts are substantially damaged. while γ R > ϕ to − ∞ and augment it, beside ∇α and ∇π as in (8), also by a strain-gradient term, referring sometimes as a conceptof non-simple materials. This may be related with an anomalous wave dispersion.In both mentioned modifications that can be even combined, ϕ is no longer component-wise quadratic. Therefore,in contrast to (8), the staggered scheme in Sect. 4 leads to non-quadratic minimization problems for which iterativesolvers have to be used. Cf. also [41] for a brief survey of various options to improve the previous model. The set of internal variables ( π , α ) from Sect. 3 is a very minimal scenario to hit basic phenomena in their simplestvariant. Actual modelling towards a geophysically relevant event should include more internal variables.Putting the model into full thermodynamical context involving temperature and heat transfer would allow toinclude, e.g., flash heating within huge earthquakes as well as the popular Dieterich-Ruina-type rate-and-state frictionmodel on the fault, cf. [37]. The energy-based ansatz (1) allows for thermodynamical consistency, i.e. the total-energyconservation, nonnegativity of temperature, and the Clausius-Duhem entropy inequality. In particular, the dissipationrate Ξ = Ξ ( q , . q ) from (3b) should be non-negative even locally for subsystems. In addition, Ξ may depend also ontemperature.Another expansion might consider a water flow in the porous rock as in [16]. It is well known that the watercontent influences both the rupture process as well as the attenuation of the seismic waves. Following Biot’s theory,the stored energy (6) is to be augmented as ϕ ( e , π , α , c ) = K E | sph e | + G E ( α ) | dev e − π | + γ DAM ( α ) + M | β ( sph e ) : I − ( c − φ ) | , (21)where c denotes the water content, and M and β are so-called Biot’s modulus and Biot’s coefficient, respectively,while φ denotes the porosity. Then the flow is governed by the Darcy law : . c = div (cid:0) m ∇ p ) with p = M (( c − φ ) − β ( sph e ) : I ) (22)where p is the pore pressure and m is the so-called hydraulic conductivity; note that p = ϕ (cid:48) c is in the position of a chemical potential and m is the mobility coefficient. Also the dissipation potential (7) can be considered c -dependent.A staggered energy-conserving time-discretisation scheme for this complex poro-thermodynamic model has beenproposed in [38]. Moreover, the porosity itself does not need to be considered fixed but may be subjected to anevolution rule similar to damage, cf. [16].Eventually, the gravitational force g in (10a) need not be considered a-priori given but rather as g = − ρ∇φ ,resulting from a gradient of the gravitational potential φ satisfying ∆φ = π g (cid:0) ρ − div ( ρ u ) (cid:1) (23)14 INETIC ENERGY DAMAGE SHEAR STORED ENERGY t = . s t = t + . s t = t + . s t = t + . s t = t + . s t = t + . s t = t + . s t = t + . s Figure 9:
Simulations of a birth of a new reverse (thrust) fault in combination with a pre-existing strike-slip faultrupturing with a certain delay, together with emission of seismic (mainly S-) waves during completion of these twosubsequenct ruptures, their propagation and creation of rather strong P-waves in the fluidic domain. The displacementagain magnified 100 × and the new mountain build during the reverse fault birth is visible. φ | | x | = ∞ = g . = . × − m kg − s − . This model of self-gravitating planet is relevant in ultra low frequency vibrations, and canbe completed with considering also centrifugal and Coriolis forces. Thus the overall bulk force is g = − ρ (cid:0) ∇φ + ω × . u + ω × ( ω × ( x + u )) (cid:1) (24)with a given angular velocity ω ∈ R and the vectorial cross product “ × ”. Some applications (and in particular those in geophysics of solid Earth) are heavily multi-scaled as far ar both timeand space concerns. Coping with this phenomenon is, numerically, very demanding and no universally applicablehints exist. Our staggered implicit scheme is well compatible with local mesh refinement which may handle spatialscales, e.g., on (and near) the fracturing surfaces (faults). Also, it is well compatible with the time-stepping refinementin (and near) the moments of ruptures and subsequent wave propagation. On the other hand, the implicit schemes ul-timately require solving large systems of algebraic equations, which brings computational restrictions. In particular,when waves are propagating through large distances, one needs a fine mesh (or high-order finite elements) every-where. For this reason, rather explicit time-discretisations are used, which ultimately, however, require the mentionedCFL condition even during periods when only slow-loading processes are on effect without any rupturing and waveemission and which makes usage of such methods limited. Here, one may think about an adaptive combination of theimplicit discretisation using rather large time-steps and admitting local mesh refinement with an explicit disretisationusing short time steps but only during relatively short events (typically minutes during earthquakes versus years orthousands of years in between them). If the 2nd-order system is implemented in a staggered way as a 1st-order system,varying of time step is easy. A rigorous combination of explicit methods with implicit discretisation of damage andMaxwellian viscoelastic rheology or poroelasticity model is not trivial and has recently be devised in [43]. The spatialscales would require still several meshes: a fine for resolving waves by explicit discretisation (possibly with specialso-called spectral finite elements [21, 25, 26]) while an only locally refined mesh for solving systems of equationsarising from the gradient-damage model.
Even solid geophysical materials (lithospheric rocks in the Earth mantle) become rather fluidic on very long timescales (millions of years). Displacements and plastic deformation can then be large even if the elastic strain issmall. A compromising model still based on small strains then involves Korteweg-like stress contribution arisingfrom the transport of internal variables, and it has been formulated in [31] and later analysed in [39]. The fullthermodynamically consistent model should be formulated at large strains, cf. [44], but its numerical implementationis very cumbersome.
A coupled model for wave source due to fracture in elastic solid (approximated by a phase-field approach) andpropagation of these waves in the layered solid and fluidic continua has been presented and tested computationallyon 2-dimensional examples. The interpretation as seismic sources by tectonic earthquakes generating seismic wavesand for their propagation was suggested. The ability of this model to capture basic phenomena occurring duringearthquakes on pre-existing faults and during the creation of new faults has been demonstrated even when usingrather basic damage-type models for the rupture.Other, more “local” and engineering, usages might be for various destructive processes in rocks or concreteconstructions adjacent to water or oil reservoirs, or in hydraulic fracture merged with natural faults (like e.g. in [24]in quasistatic variant).The algorithmically efficient staggered time discretisation scheme has been used. This was combined with thesimplest P1 finite elements for the space discretisation. Yet, it is well known that this discretisation is not optimal forelasticity and more sophisticated spectral elements [21–23, 25, 26] or discontinous P2-elements [45] are implementedin this context. 16ventually, various expansions of this basic model have been outlined. In fact, Sections 2–4 already introduceanother internal variable (the inelastic strain) that should be considered for long-time processes as creep or as healingduring long lasting periods in between earthquakes, although it was not implemented in the presented computationalsimulations in Sect. 5. R EFERENCES [1] L. Ambrosio and V. M. Tortorelli. Approximation of free discontinuity problems.
Boll. Unione Mat. Italiana ,6-B:105–123, 1992.[2] Z.P. Baˇzant and M. Jir´asek. Softening-induced dynamic localization instability: seismic damage in frames.
J.Engr. Mech. , 122:1149–1158, 1996.[3] A. Bedford.
Hamilton’s Principle in Continuum Mechanics . Pitman, Boston, 1985.[4] Y. Ben-Zion. Properties of seismic fault zone waves and their utility for imaging low-velocity structures.
J.Geophys. Res. , 103:12,567–12,585, 1998.[5] Y. Ben-Zion. Dynamic ruptures in recent models of earthquake faults.
J. Mech. Phys. Solids , 49:2209–2244,2001.[6] Y. Ben-Zion and J.-P. Ampuero. Seismic radiation from regions sustaining material damage.
Geophys. J. Int. ,178:1351–1356, 2009.[7] D. V. Boger. A highly elastic constant-viscosity fluid.
J. Non-Newtonian Fluid Mechanics , 3:87–91, 1977.[8] B. Bourdin, C.J. Larsen, and C.L. Richardson. A time-discrete model for dynamic fracture based on crackregularization.
Int. J. of Fracture , 10:133–143, 2011.[9] B. Bourdin, J.-J. Marigo, C. Maurini, and P. Sicsic. Morphogenesis and propagation of complex cracks inducedby thermal shocks.
Phys. Rev. Lett. , 112:014301, 2014.[10] E. Choi, E. Tan, L.L. Lavier, and V.M. Calo. DynEarthSol2D: An efficient unstructured finite element methodto study long-term tectonic deformation.
J. Geophys. Res. Solid Earth , 118:2429–2444, 2013.[11] R. Courant, K. Friedrichs, and H. Lewy. ¨Uber die partiellen Differenzengleichungen der mathematischen Physik.
Math. Annalen , 100:32–74, 1928.[12] M. Focardi. On the variational approximation of free-discontinuity problems in the vectorial case.
Math. ModelsMethods Appl. Sci. , 11:663–684, 2001.[13] A. Green and P. Naghdi. A general theory of an elastic-plastic continuum.
Arch. Rational Mech. Anal. , 18:251–281, 1965.[14] A.A. Griffith. The phenomena of rupture and flow in solids.
Phil. Trans. Royal Soc. London, A , 221:163–198,1921.[15] B. Halphen and Q.S. Nguyen. Sur les mat´eriaux standards g´en´eralis´es.
J. M´ecanique , 14:39–63, 1975.[16] Y. Hamiel, V. Lyakhovsky, and A. Agnon. Coupled evolution of damage and porosity in poroelastic media:theory and applications to deformation of porous rocks.
Geophys. J. Int. , 156:701–713, 2004.[17] W.R. Hamilton. On a general method in dynamics, part II.
Phil. Trans. Royal Soc. , pages 247–308, 1834.[18] R. A. Harris et al. The SCEC/USGS dynamic earthquake rupture code verification exercise.
Seismological Res.Lett. , 80:119–126, 2009.[19] M. Hofacker and C. Miehe. Continuum phase field modeling of dynamic fracture: variational principles andstaggered FE implementation.
Intl. J. Fracture , 178:113–129, 2012.1720] Y. Huang, J.-P. Ampuero, and D. V. Helmberger. Earthquake ruptures modulated by waves in damaged faultzones.
J. of Geophysical Research: Solid Earth , B9:3133–3154, 2014.[21] Y. Kaneko, N. Lapusta, and J.-P. Ampuero. Spectral element modeling of spontaneous earthquake rupture on rateand state faults: Effect of velocity-strengthening friction at shallow depths.
J. Geophysical Res. , 113:B09317,2008.[22] M. K¨aser and M. Dumbser. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstruc-tured meshes I: The two-dimensional isotropic case with external source terms.
Geophys. J. Int. , 166:855–877,2006.[23] M. K¨aser, M. Dumbser, J. de la Puente, and H. Igel. An arbitrary high-order discontinuous Galerkin method forelastic waves on unstructured meshes III: Viscoelastic attenuation.
Geophys. J. Int. , 168:224–242, 2007.[24] A.R. Khoei, M. Vahab, and M. Hirmand. Modeling the interaction between fluid-driven fracture and naturalfault using an enriched-FEM technique.
Intl. J. Fracture , 197:1–24, 2016.[25] D. Komatitsch and J. Tromp. Spectral-element simulations of global seismic wave propagation - I. validation.
Geophys. J. Int. , 149:390–412, 2002.[26] D. Komatitsch and J. Tromp. Spectral-element simulations of global seismic wave propagation - II. three-dimensional models, oceans, rotation and self-gravitation.
Geophys. J. Int. , 150:303–318, 2002.[27] M. Kruˇz´ık and T. Roub´ıˇcek.
Mathematical Methods in Continuum Mechanics of Solids . Springer, Switzerland,2019.[28] C.J. Larsen, C. Ortner, and E. S¨uli. Existence of solution to a regularized model of dynamic fracture.
Math.Models Meth. Appl. Sci. , 20:1021–1048, 2010.[29] T. Lay and T. C. Wallace.
Modern global seismology . Acad. Press, San Diego, 1995.[30] V. Lyakhovsky and Y. Ben-Zion. A Continuum Damage-Breakage Faulting Model and Solid-Granular Transi-tions.
Pure Appl. Geophysics , 171:3099–3123, 2014.[31] V. Lyakhovsky and Y. Ben-Zion. Damage-breakage rheology model and solid-granular transition near brittleinstability.
J. Mech. Phys. Solids , 64:184–197, 2014.[32] V. Lyakhovsky, Y. Hamiel, J.-P. Ampuero, and Y. Ben-Zion. Non-linear damage rheology and wave resonancein rocks.
Geophys. J. Int. , 178:910–920, 2009.[33] V. Lyakhovsky, Y. Hamiel, and Y. Ben-Zion. A non-local visco-elastic damage model and dynamic fracturing.
J. Mech. Phys. Solids , 59:1752–1776, 2011.[34] V. Lyakhovsky and V.P. Myasnikov. On the behavior of elastic cracked solid.
Phys. Solid Earth , 10:71–75,1984.[35] A. Mielke and T. Roub´ıˇcek.
Rate-Independent Systems – Theory and Application . Springer, New York, 2015.[36] C. Pelties, J. de la Puente, J.-P. Ampuero, G. B. Brietzke, and M. K¨aser. Three-dimensional dynamic rupturesimulation with a high-order discontinuous galerkin method on unstructured tetrahedral meshes.
J. Geophys.Res. , 117:B02309, 2012.[37] T. Roub´ıˇcek. A note about the rate-and-state-dependent friction model in a thermodynamical framework of theBiot-type equation.
Geophys. J. Intl. , 199:286–295, 2014.[38] T. Roub´ıˇcek. An energy-conserving time-discretisation scheme for poroelastic media with phase-field fractureemitting waves and heat.
Disc. Cont. Dynam. Syst. S , 10:867–893, 2017.[39] T. Roub´ıˇcek. Geophysical models of heat and fluid flow in damageable poro-elastic continua.
Cont. Mech.Thermodyn. , 29:625–646, 2017. 1840] T. Roub´ıˇcek. Seismic waves and earthquakes in a global monolithic model.
Cont. Mech. Thermodynam. ,30:709–729, 2018.[41] T. Roub´ıˇcek. Models of dynamic damage and phase-field fracture, and their various time discretisations. In J.-F.Rodrigues and M. Hintterm¨uller., editors,
Topics in Applied Analysis and Optimisation , CIM Series in Math.Sci. Springer, in print. Preprint arXiv 1906.04110.[42] T. Roub´ıˇcek and C. G. Panagiotopoulos. Energy-conserving time-discretisation of abstract dynamical problemswith applications in continuum mechanics of solids.
Numer. Funct. Anal. Optim. , 38:1143–1172, 2017.[43] T. Roub´ıˇcek, C.G. Panagiotopoulos, and C. Tsogka. Explicit time-discretisation of elastodynamics with someinelastic processes at small strains . 2019. Preprint arXiv no.1903.11654..[44] T. Roub´ıˇcek and U. Stefanelli. Thermodynamics of elastoplastic porous rocks at large strains towards earthquakemodeling.
SIAM J. Appl. Math. , 78:2597–2625, 2018.[45] J. Tago, V. M. Cruz-Atienza, J. Virieux, V. Etienne, and F. J. Snchez-Sesma. A 3D hp-adaptive discontinuousGalerkin method for modeling earthquake dynamics.
J. Geophys. Res. , 117:B09312, 2012.[46] N. Tosi, O. ˇCadek, and Z. Martinec. Subducted slabs and lateral viscosity variations: effects on the long-wavelength geoid.