Seismic waves and earthquakes in a global monolithic model
Tomáš Roubíček







Seismic waves and earthquakes in
a global monolithic model.






Tomáš Roubíček^1,2






^1 Mathematical Institute, Charles University,
Sokolovská 83, CZ-18675 Praha 8, Czech Republic

^2 
Institute of Thermomechanics, CAS,


Dolejškova 5, 
CZ-18200 Praha 8, Czech Republic 






Abstract:
The philosophy that a single “monolithic” model can “asymptotically” 
replace and couple in a simple elegant way several specialized models relevant 
on various Earth layers is presented and, in special situations, also 
rigorously justified. In particular, global seismicity and tectonics is 
coupled to capture e.g. (here by a simplified model) ruptures of 
lithospheric faults generating seismic waves which then propagate through the 
solid-like mantle and inner core both as shear (S) or pressure (P) waves, 
while S-waves are suppressed in the fluidic outer core and also in the oceans. 
The “monolithic-type” models have the capacity to describe all the 
mentioned features globally in a unified way together with corresponding 
interfacial conditions implicitly involved, only when scaling its parameters 
appropriately in different Earth's layers. Coupling of seismic waves with 
seismic sources due to tectonic events is thus an automatic side effect. 
The global ansatz is here based, rather for an illustration, only on a 
relatively simple Jeffreys' viscoelastic damageable material at small strains 
whose various scaling (limits) can lead to Boger's viscoelastic fluid or even 
to purely elastic (inviscid) fluid. Self-induced gravity field, Coriolis, 
centrifugal, and tidal forces are counted in our global model, as well. 
The rigorous mathematical analysis as far as the existence of
solutions, convergence of the mentioned scalings, and energy conservation 
is briefly presented.



AMS Classification: 
35K51, 

35L20, 






35Q86, 
74J10, 
74R20, 
74F10, 


76N17, 




86A15, 
86A17. 



Key words: waves; global seismicity; tectonic earthquakes; 
mathematical models; energy conservation; scaling; 
convergence proofs, existence of weak solutions. 

 





§ INTRODUCTION



Global geophysical models typically (have to) deal with several very 
different phenomena and couple various models due to the layered character 
of our planet Earth as well as our Moon and many other terrestrial planets 
or other moons. This paper wants to demonstrate (and, in special
situations, also rigorously justify) the philosophy that a single 
“monolithic” model can replace several specialized models coupled together.
Such a single model can be  easier to implement on computers algorithmically. 
Of course, computationally, such a monolithic-type model may be not 
always easier to produce really relevant simulations at computers 
we have at our disposal nowadays. As purely seismic global 3D models 
are already treatable as well as their coupling with earthquake source 
at least locally, cf. e.g. 
<cit.>
as well as <cit.>, 
respectively, there 
is a hope that the coupling monolithic approach may become more amenable 
in future within ever increasing computer efficiency.


The phenomena we have in mind in this paper involve global seismicity and 
tectonics. In particular the latter involves e.g. ruptures of 
lithospheric faults generating seismic waves which then propagate through the 
solid-like silicate mantle and iron-nickel inner core both in shear (S) or 
pressure (P) modes. In contrast to the P-waves (also called primary or
compressional waves), the S-waves (also called secondary waves) are suppressed 
in the fluidic iron-nickel outer core and also in the water oceans 
(where P-waves emitted from the earthquakes in the crust may manifest as 
Tsunami at the end). 

Typically, Maxwell-type rheologies are used in geophysical models
of solid 
mantle to capture long-term creep effects up to 10^5 yrs. Sometimes, 
more attenuation of Kelvin-Voigt type is also involved, which leads 
to the Jeffreys rheology <cit.>, cf. also <cit.>.
This seems more realistic in particular because it covers (in the limit)
also the Kelvin-Voigt model applied to the volumetric 
strain whereas the pure Maxwell rheology allowing for big creep
during long geological periods is not a relevant effect in the volumetric part.
We therefore take the Jeffreys model as a basic global “monolithic” 
ansatz and, in various limits in the deviatoric and the volumetric parts, 
we model different parts of the planet Earth.

Respecting the solid parts of the model, we use the 
Lagrangian description, i.e. here all equations are formulated 
terms of displacements rather than velocities, while the reference
and the actual space configurations automatically coincide with each other in our
small strain (and small displacement) ansatz, which is well relevant
in geophysical short-time scales of seismic events. 

In the solid-like part, various inelastic processes are considered 
to model tectonic earthquakes on lithospheric faults together with
long healing periods in between them, as well as aseismic slips, and
various other phenomena. To this goal, many internal variables
may be involved as aging/damage, inelastic strain, porosity, water content, 
breakage, and temperature, cf. <cit.>.
On the other hand, those sophisticated models are focused on 
rather local events around the tectonic faults without ambitions 
to be directly coupled with the global seismicity. 
Cf. also <cit.> for models of 
rupturing lithospheric faults and, in particular, a relation with 
the popular Dieterich-Ruina rate-and-state friction model. Here, rather for 
the lucidity of the exposition, we reduce the set of internal variables 
to only one scalar variable, namely damage/aging, which 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 overall global model. Also, this simple model already 
will well illustrate mathematical difficulties related to nonlinearities 
in the solid parts coupled with linear but possibly hyperbolic fluidic regions.

Let us emphasize that usual models are focused only either on propagation of 
seismic waves along the whole globe while their source is considered given, or 
on description of seismic sources due to tectonic events, 
but not their mutual coupling.

If a coupling is considered, then is concerns rather local models 
not considering the layered structure of the whole planet, cf. e.g.  
<cit.>.
The reality ultimately 
captures very different mechanical properties of different 
layers of the Earth, in particular the mantle and the inner core which are 
solid from the short-time scales versus the outer core and the oceans 
which are fluidic even on the short-time scales.

 
The goals of this article are:



  )to propose a model that might
capture simultaneously propagation of seismic waves over the whole planet
and their nonlinearly behaved sources (like ruptures of tectonic faults, 
here modelled only in a very simplified way for a relatively lucid 
illustration of the model procedure), both mutually coupled. 


  )
by proper scaling to approximate viscoelastic Boger-type <cit.>
fluids that is relevant in outer core and in the oceans (with a very low 
viscosity) where S-waves can then only slightly penetrate the outer core or
the oceans but are fast attenuated, while P-waves are only refracted.


  )
by limiting further the viscosity in the outer core or the oceans to zero, 
further to approximate this viscoelastic fluid towards  elastic (completely
inviscid) fluid, 

respecting the phenomenon that 
S-waves cannot penetrate into these fluidic regions and are fully reflected 
on the interfaces between the outer core and the mantle 
(=Gutenberg's discontinuity) and the inner core
or on the ocean beds, i.e. between  and , while P-waves 
propagate through these interfaces, being both refracted and reflected on them. 


  )
perform the rigorous analysis as far as the existence of solutions,
a-priori estimates in specific norms, and convergence towards other models
that justifies the particular models, their energetics and asymptotics,
and can support numerical stability and convergence when discretised
and implemented on computers.  


The plan of this article is the following: In Section <ref>, the general 
monolithic model is introduced and also specialized for viscoelastic 
linear isotropic material undergoing possibly damage in the elastic shear part.
In Sections <ref> and <ref>, we formulate the limit toward 
viscoelastic and purely elastic fluids, respectively. This limit concerns
only the sudbomain , i.e. the outer core and oceans. To make the 
presentation more lucid and readable also for the geophysical specialist 
outside mathematical-analyst community, the mathematical analysis of the 
models and necessarily a bit technical proofs of the statements formulated in 
Sections <ref>–<ref> are intentionally presented as Appendix later in 
Section <ref>. For it, we use a conceptually constructive approximation
by the Galerkin method, which is in some variant also used in geophysical 
literature, cf. e.g. <cit.>.





§ PHILOSOPHY AND AN EXAMPLE OF A MONOLITHIC MODEL



π

Rather for illustration of our main focus to coupling of solid and fluidic
models, we ignore most of the above mentioned internal variables and 
keep only the scalar-value aging/damage variable, denoted by α,
and then we consider the mentioned Jeffreys rheology combined with damage 
influencing only the elastic but not viscous part. 
The differential relation governing this rheology is of the type 
“σ+σ= e+ e” with some specific 
coefficients (here ignored), some of them being possibly subject to damage. 
Here and in what follows, the dot-notation (·)^
stands for the time derivative ∂/∂ t.
Alternatively, to make the damage dependence more 
lucid, one can use system of 
first-order equations in time when implementing the concept of
internal variables, i.e. here introducing a Maxwellian-type 
creep strain  

and making a Green-Naghdi <cit.>
additive decomposition of the total strain 

    e(u)=e_ el+.

So, considering also the damage/aging variable, altogether our set of 
internal variables will be (,α).

The main ingredients of the model are then, beside the mass density 
ϱ=ϱ(x), the specific stored energy and the dissipation 
potential as a property of the material 


    φ=φ(x,u,e_ el,α,∇α,ϕ,∇ϕ),
    
       ζ=ζ(x, e_ el,,α).

Note that φ and ζ naturally do not explicitly 
depend on the total strain e(u) but 
rather on the elastic strain e_ el and possibly on the 
creep strain  and their rates; in fact, φ is independent
of  because naturally no hardening-like effects are relevant in
geophysical models. In what follows, we will often omit 
the explicit dependence on x variable.

For readers' convenience, Table 1 summarizes the main nomenclature used 
in this paper.





To formulate the equations for the dynamics of displacement u, it is more 
natural in the Kelvin-Voigt model (which is a part
of the Jeffreys rheology) to express the stored-energy and dissipation 
potentials 
rather in terms of the total strain e=e(u) and the creep  from 
e-decompose as


    φ(x,u,e,,α,∇α,ϕ,∇ϕ)
    =φ(x,u,e-,α,∇α,ϕ,∇ϕ),
    
       ζ( e,,α
    )
    =
    ζ( e-,,α
    ).

In terms of these potentials, our system in the abstract form 
is: 

    ϱ u+φ_u'(∇ϕ)
    - div (
    ζ_ e'(e( u),)+φ_e'(e(u),,α))
    
       =f__ COR( u)
    
    
          in  Ω,      
    
    
    ζ_'(e( u),)+φ_'(e(u),,α)
    
       =0
    
          in  Ω,
    
    
    ζ_α'(α)
    +φ_α'(e(u),,α)
    
    - div(
    
    φ_∇α'(∇α))   =0      in  Ω,
    
    
    
    φ_ϕ'(u,ϕ)- div(φ_∇ϕ'(∇ϕ))
       =ϱ_ ext(t)
          in  ^3,

where (·)' denotes the (partial) derivatives for which 
we already have assumed a certain particular form of phi-zeta
so that the partial derivatives of φ and ζ do not
depend on the full list of arguments, thus e.g. instead of 
φ_e'(u,e(u),,α,∇α,ϕ,∇ϕ) we wrote only 
φ_e'(e(u),,α), etc. Note that, in system-u, we 
used also a bulk force f__ COR( u) which may not come from the
dissipation potential by a classical way
f__ COR( u)=ζ_ u, which in particular
the case of the Coriolis (pseudo) force which itself does not 
vanish but nevertheless the dissipation rate f__ COR( u) u=0.

Rheological response under the volume variation
and under the shear may (and do) substantially differ from each other.
To distinguish these geometrical aspects,  the total strain is decomposed to
the spherical (also called hydrostatic or volumetric) 
strainstrain!spherical 

and the deviatoric (also called shear) 
strain:strain!deviatoric 

    e= sph e+ dev e    with     sph e= tr e/3    and     dev e=e- tr e/3,

where  is the identity matrix and “tr” denotes the trace of
a matrix. Note that the deviatoric and the spherical strains 
from decomposition-of-strain are orthogonal to each other.
In terms of this decomposition, for example the isotropic elastic 
(Lamé) material at small strains has the quadratic stored energy

    φ__(e)=1/2λ( tr e)^2+G__ E|e|^2=
    1/2(λ+2/3G__ E)( tr e)^2
    +G__ E|e- tr e/3|^2
    
               =
    1/2(3λ+2G__ E)| sph e|^2
    +G__ E| dev e|^2
    
               =3/2K__ E| sph e|^2+G__ E| dev e|^2    with  K__ E=λ+2/3G__ E.

where we used | sph e|^2=( tr e/3)^2:=
( tr e/3)^23=( tr e)^2/3 and 
the mentioned orthogonality of sph e and dev e,

where 
K__ E is called the bulk modulus and λ is 
the (first) Lamé coefficient while G__ E is the shear modulus
(= the second Lamé coefficient).
Then the corresponding stress writes as:

    φ__'(e)   =λ( tr e)
    
    +2G__ E e=:σ_ dev+σ_ sph
          with   
    σ_ dev=2G__ E dev e    and   σ_ sph=3K__ E sph e.


We will use this decomposition both for elastic and for the viscous parts
of our model and let the elastic bulk modulus G__ E and thus also 
φ dependent on damage. Yet, this dependence would bring mathematical 
difficulties (like lack of the integrability of the term α G__E(α)
| dev e_ el|^2 that would arise in system-alpha-w-, for 
example). Here, various modifications of the Lamé-type model can help. 
One option is to consider higher strain gradients, i.e. the concept
of nonsimple continua, cf. also <cit.>. Other (conceptually
and technically simpler) option is to modify the stored energy by considering
an non-quadratic e↦__ E(α, dev e) instead
of e↦ G__E(α)| dev e|^2. A canonical option one may have
in mind is

    __ E
    (α,e_ el)=__ E
    (α, dev e_ el)=
    G__ E(α)| dev e_ el|^2/√(1+| dev e_ el|^2/E__ M^2)


for (presumably large) regularizing parameter E__ M. 
In accord with Figure <ref>, in KV-limit+ and
in what follows, we quite naturally assume that __ E does not depend 
on the spherical part of the strain which itself usually cannot make any
damage in the rock. The meaning of E__ M in (<ref>)
is clear from realizing that, if the deviatoric 
strain is substantially smaller in the sense 
| dev e_ el|<<E__ M, then the difference from the original 
Lameé-type model is negligible, i.e. the corresponding stress contribution 
∂_e 

__ E(α, e_ el)∼ 2G__ E(α)
 dev e_ el.
On the other hand, this modification is effective for large strains and makes 
e_ el↦φ(x,u,e_ el,α,∇α,ϕ,∇ϕ) 
nonquadratic but still keeps convex now with at most linear growth,
so that the elastic part of deviatoric stress is always below 
G__ E(α)E__ M. Thus, for E__ M large, this is effectively
not a substantial (and physically well acceptable) modification 
of the original model. 





Thus we consider the material stored energy and the 
material dissipation-force potential from phi-zeta governing the 
problem as 


    φ=φ(x,u,e_ el,α,∇α,ϕ,∇ϕ)=
    3/2K__ E| sph e_ el|^2+
    
    
    
    __ E(α, dev e_ el)
    
    
       
    
                          +κ/2|∇α|^2+ϱϕ
    +ϱ∇ϕ·u+1/8 g|∇ϕ|^2   
    
                          +1/2ϱ(|ω·(x+u)|^2-|ω|^2|x+u|^2)
       if x∈Ω,
    
    =
    φ(∇ϕ)=1/8 g|∇ϕ|^2   if x∈^3∖Ω,
    
    
       ζ=ζ(x, e_ el,,α
    )=3/2K__ KV| sph  e_ el|^2
    +G__ KV| dev  e_ el|^2
    
                               
    +3/2K__ MX| sph |^2
    +G__ MX| dev |^2
    +η(α)

with g≐6.674×10^-11m^3kg^-1s^-2 the gravitational constant. 
For the ϕ-terms in def-of-phi see e.g. 
<cit.>. 
All the coefficients K's, G's, and κ's 
are naturally considered defined on Ω and x-dependent (which is not 
explicitly written in def-of-phi-zeta for brevity).
Also the mass density is considered  x-dependent
and, because of system4, defined on the whole Universe ^3 
by putting ϱ(x)=0 if x∈^3∖Ω.
Keeping K__ E constant in def-of-phi reflects the phenomenon that
compression does not cause damage and the pure tension is not 
a relevant mode in geophysical models, so only shear does cause
damage because of dependence of the elastic shear modulus 
G__ E in KV-limit+ and def-of-phi on α.
The last term in def-of-phi is the potential of the centrifugal
force ϱω×(ω×(x+u))
and, noteworthy, it violates coercivity due to the
term -1/2ϱ|ω|^2|x+u|^2, which reflects the real phenomenon
that, given the angular velocity ω, the centrifugal force
can indeed inflate the body in an unlimited way in the direction orthogonal
to the rotation axis. Of course, in reality, this is either not 
observed in planetary/satellite bodies or the angular velocity cannot be taken
apriori fixed in some large-type bodies. 


For shorter notation, we define the 4th-order tensors corresponding
to the isotropic material expressed in terms of the (K,G)-moduli in 
def-of-phi-zeta, namely


    [(α,e)]_ijkl= K__ Eδ_ijδ_kl
    +∂_ee^__ E(α,e)(δ_ikδ_jl+δ_ilδ_jk-2δ_ijδ_kl/3),
    
       [__ MX]_ijkl=K__ MXδ_ijδ_kl
    +G__ MX(δ_ikδ_jl+δ_ilδ_jk-2δ_ijδ_kl/3),
    
       [__ KV]_ijkl=K__ KVδ_ijδ_kl
    +G__ KV(δ_ikδ_jl+δ_ilδ_jk-2δ_ijδ_kl/3),

where δ's denote the Kronecker symbol.
In what follows, we will take specifically 
f__ COR( u)=2ϱω× u
which is the standard Coriolis force with
ω the given vector (angular velocity
related with the constant rotation of the Earth with respect
to the inertial system).
The general system system then takes the more specific form:  

    ϱ u- div σ
    
    +f=0   with  σ=__ KV e_ el
    +(α, dev e_ el)e_ el 
       and  
    f=ϱ(ω×(ω×(x+u))+2ω× u+∇ϕ)      in  Ω,
    
       
    __ MX
    
    =σ
          in  Ω,
    
       η(α)
    
    +
    
    ∂_α__ E(α, dev e_ el)
    - div(κ∇α)∋0
    
          in  Ω,
    
       
     div(1/4 g∇ϕ+ϱ u)=ϱ
    +ϱ_ ext(t)
          in  ^3.

The bulk force f in system-u involves the centrifugal force, the 
Coriolis force, and the 
gravity force due to the self-induced gravity field. This last force plays 
a certain role in ultra-low frequency (i.e. very long wavelength) seismic 
waves, cf. e.g. <cit.>. 
For such right-hand side in system-u+ together with 
system4 see e.g. <cit.>. Actually,
this self-gravity interaction results after a certain linearization 
of the system originally written at large strains, cf. also 
<cit.>. The external time-varying 
mass density ϱ_ ext^=ϱ_ ext^(x,t) occurring in 
system4 allows for involvement of tidal forces arising from 
other bodies (stars, planets, moons) moving around.
The notation η in system-alpha stands for a subdifferential
of the convex function η which standardly generalizes the usual
derivative if η is non-smooth and, by definition, the inclusion 
“∋” there actually means, for any v and a.a. (t,x)∈ Q, the inequality 

    (v-α(t,x))
    (∂_α__ E(x,α(t,x), dev e_ el(t,x))
    
    - div(κ∇α(t,x)))
    +η(v)≥η(α(t,x)).





To facilitate spontaneous rupture on some weak narrow layers (pre-existing 
faults in the solid crust) that leads to earthquakes, some 
“enough pronounced” nonconvexity of φ is necessary. One option was 
devised in <cit.>, introducing non-convexity in terms
of e_ el if damage develops enough, and used in a series of papers, 
see e.g. <cit.> and references therein. 
Mathematically rigorous setting of this model requires however 
some higher-order strain gradients (the so-called 2nd-grade nonsimple
materials) and is rather technical, cf. <cit.>. 
Alternative nonconvexity can be implemented through 
__ E(·,e) causing 
possible weakening when damage develops, being opposite to so-called cohesive 
damage. To avoid the mentioned mathematical technicalities, we have
chosen the latter option.


To see the energetics behind the system system, we test the particular 
equations in system by u, π, α, and ϕ, 
respectively. 

Thus we obtain (at least formally) the energy balance 

    /ṭ(∫_Ωϱ/2| u|^2 x̣+
    ∫_^3φ(u,e(u),,α,∇α,ϕ,∇ϕ) x̣)      
    +∫_Ωξ(e( u),,α) x̣
    =∫_^3ϱ_ ext(t)ϕ x̣ ,

where ξ=ξ( e,,α) is the dissipation rate related 
with the dissipation potential ζ=ζ( e,,α) by 
ξ( e,,α)
=( e,,α)^⊤∂ξ( e,,α).


Using def-of-phi-zeta and C-D, we can make it 
more specific for the system system+ when using also the calculus

    ∫_Ωϱ∇ϕ u- div(ϱ u)ϕ x̣
    =∫_Ωϱ∇ϕ u+ϱ u∇ϕ x̣
    =/ṭ∫_Ωϱ∇ϕ u x̣.
 
and realizing that, in view of  def-of-zeta- and def-of-zeta,

    ξ( e,,α)   =3K__ KV| sph( e-)|^2
    +2G__ KV| dev( e-)|^2
    
       
      +3K__ MX| sph |^2
    +2G__ MX| dev |^2
    +α η'(α)
    
       
    =__ KV( e-π)∷( e-π)
    +__ MXπ∷π+αη'(α)

when using also the notation (<ref>b,c).
Thus engr-abstract reads more speficially here as 

    /ṭ(∫_Ωϱ/2| u|^2+3/2K__ E| sph(e(u)-π)|^2
    +
    
    __ E(α, dev(e(u)-π))
    +ϱϕ+ϱ∇ϕ u
    
       
    +κ/2|∇α|^2
    +1/2ϱ(|ω·(x+u)|^2-|ω|^2|x+u|^2)
     x̣
    +∫_^31/8 g|∇ϕ|^2 x̣)
    
    
       
    +∫_Ω__ KV(e( u)-π)∷ (e( u)-π)
    +__ MXπ∷π+α∂η(α) x̣
    =∫_^3ϱ_ ext(t)ϕ x̣.

Here one should naturally assume that η is (naturally) nonsmooth possibly 
only at α=0 so that the dissipation rate 
α∂η(α) is actually 
a single-valued function.
The Coriolis force does not occur in engr
because ω× u is always orthogonal to 
u so that (ω× u) u=0. This is a well 
known effect that the Coriolis (pseudo) force does not make any work.

Denoting n⃗ the (unit) normal to the boundary ∂Ω of 
Ω, i.e. the Earth surface, we complete the system by natural initial 
and boundary conditions: 

    u|_t=0=u_0,     u|_t=0=v_0,    π|_t=0=π_0,    α|_t=0=α_0.
    
       σn⃗=0         on ∂Ω,
        ϕ(∞)=0,   and  
       
    (κ∇α)·n⃗=0  on ∂Ω,

Actually, it would make a good sense to consider also some period instead of initial conditions,
but the analysis of such a problem would be more difficult.

It is a conventional modelling property that α ranges the interval 
[0,1], with α=0 meaning no damage while α=1 corresponds to a 
completely disintegrated rock. To make the damage model relatively simple 
without causing unnecessary analytic complications (leading e.g. to 
implementing the concept of nonsimple materials), the simplest modelling 
assumption that ensures the mentioned constraint 0≤α≤1 is 

    ∂_α__ E(0,e)=0=∂_α__ E(1,e).
  
This facilitates the analysis of the model, sketched in Appendix below. Without 
going into details here, we can here summarize the main theoretical result as:


The initial-boundary-value problem system+–BC-IC admits a weak 
solution in the sense of Definition <ref> below.
This solution conserves energy in the sense that engr holds when
integrated over any time interval [0,t] with 0<t≤ T.



Some other modification of the “solid” model used in geophysical
literature modifies the evolution of the creep  to be 
activated in the spirit what it standardly used in plasticity. Then
 ζ should be augmented by a term like | dev| and,
to facilitate mathematical analysis, then φ is to be augmented  
by the term |∇ dev|^2. Then  is called an inelastic
strain, rather than a Maxwellian creep. This model may be relevant in particular



to provide an additional dissipation of
energy important during big earthquakes and damage-dependent 
yield stress that may facilitate fast rupture.



In the standard presentation of the selfgraviting model,
the evolving potential ϕ is rather a perturbation of 
another gravitational potential (constant in time), resulted from 
a steady-state equilibrium configuration at a chosen initial time
with the Coriolis force 

zero and the centrifugal force 
ϱω×(ω×x). The sum of these two is 

sometimes referred as a geopotential. Then another force occurs in 
system-u+, causing the body Ω to be pre-stressed. Cf. e.g. 
<cit.>. 










§ TOWARDS BOGER VISCOELASTIC FLUID IN THE OUTER CORE AND OCEANS




A standard specification of so-called Boger's viscoelastic fluids 
is as constant-viscosity elastic (non-Newtonian) fluids that 
behave as both liquids and solids. Originally, this concept was devised rather 
for materials like dilute polymer solutions. 
Yet, such (idealized) material may model the fluid outer core in the Earth,
which is a 2200 km thick layer between the (rather solid) inner 
core of the radius about 1300 km and the (rather solid) mantle, 
cf. Fig. <ref>. The elasticity in the volumetric part is essential 
because it allows for a propagation of P-waves, which is a well documented 
phenomenon in the outer core, while S-waves are practically not penetrating 
this region because of its fluidic character as far as the shear response 
concerns, cf. Figure <ref>.

Knowing quite reliably the speed v__ P of the P-waves 
and the mass density ϱ depending on the depth, 
the elastic bulk modulus K__ E=K__ E(x) in the outer core is easy to 
be seen from the formula v__ P=√(M/ϱ) where 
M=λ+2G__ E=K__ E+4/3G__ E is the so-called P-wave modulus;
here G__ E=G__ E(α) in KV-limit+. In the 
elastic fluid G__ E=0 so that M=K__ E. One can indeed determine simply as
K__ E=ϱ v__ P^2 when knowing ϱ and v__ P. 
According the reliably documented speed of the seismic P-waves
varying from 8 km/s (top) to 10 km/s (bottom) and the mass density 
varying from 10000 kg/m^3 (top) to 12000 kg/m^3 (bottom),
the elastic bulk modulus ranges from K__ E=640GPa (top) to 
1200 GPa (bottom).

For oceanic layers, it is relevant that the bulk modulus of seawater is about 
K__ E=M=2.3GPa (increasing with pressure) which is, e.g., about 
80× smaller than in steel and about 5× smaller than in the crust, 
and thus water is remarkably elastically compressible. 
Modelling oceans as elastic fluid is thus relevant for propagation
seismic P-waves, which may manifest themselves on the surface by the 
Tsunami waves, cf. e.g. 
<cit.>. 
Note that, counting  the speed of P-waves close to the 
sea surface where ϱ∼10^3 kg/m^3 is 
v__ P=√(M/ϱ)∼1.5km/s, which is the sound speed,
so underwater acoustics exploits the same mechanism as seismic P-waves.



In addition to the elastic response, there is a (rather small) viscosity 
(K__ KV,G__ KV) 
in the fluidic outer core and in oceans, too. In the outer core, some data 
indicates it of the order 10^-2Pa s, cf. <cit.>, while some 
uncertainty and big variation within depth 10^4Pa s
seems also documented in literature, 
see e.g. <cit.>. The viscosity of the water
is even lower, of the order 10^-3Pa s, varying with temperature and 
pressure. 

Here we can come to a viscoelastic (Boger's) fluid as a limit from the 
previous damageable Jeffreys' rheology when sending G__ E(α)→0 
so that, in particular, the resulting material naturally will 
not be subjected to any damage, and simultaneously 
K__ MX→∞, and G__ MX→∞.
As a result, π→0 and eventually system-pi is to hold only
on  instead of the whole Ω while
system-alpha still holds on the whole Ω
but 

_α__ E(α,e_ el)
is zero on 
and values of the damage α are actually irrelevant 
 system-u+ on .

More specifically, the system that results from system+
by this limit procedure looks as:

    ϱ u- div(σ_ sph+σ_ dev)
    
    
    +ϱ(ω×(ω×(x+u))+2ω× u+∇ϕ)=0      in  Ω
        with  σ_ sph=3K__ KV sph  e_ el
    +3K__ E sph
     e_ el        in  ,
    
    3K__ KV sph e( u)
    +3K__ E sph e(u)   in  ,
        and   σ_ dev=
    2G__ KV dev  e_ el
    +
    
    ∂_e__ E(α, dev e_ el)
    
       in  ,
    
    2G__ KV sph e( u)   in  ,
    
       
    __ MX
    
    =σ_ sph+σ_ dev
          in  ,
    
    
       η(α)
    
    +
    
    
    ∂_α G__ E(α, dev e_ el)
    
    - div(κ∇α)∋0
    
          in  ,
    
    
    
    
    
       
     div(1/4 g∇ϕ+ϱ u)
    =ϱ+ϱ_ ext(t)
    
          in  ^3.

The interface conditions for displacement/stress on the interior boundaries 
between the mantle and the 
outer core and the inner and the outer core as well between the 
mantle (crust) and the oceans are automatically involved as system-u++
holds on the whole Ω. In particular, the stress-vector equilibrium and 
continuity of the displacement across these interfaces are automatically
involved and does not need to be written explicitly. On the other hand,
note in particular that both π and α in system++
are now needed and defined only in the solid part .








Therefore, the condition for the “flux” of α is now to be prescribed 
on all interior boundaries (i.e. on the mantle/core and mantle/oceans and
inner/outer-core interfaces). More specifically,
BC-IC-alpha is to be replaced by

    (κ∇α)·n⃗=0  on ∂.


We can justify this limit-passage scenario rigorously. To this goal, let us 
choose 

    G__ MX,(x)=G__ MX(x)/,   
     G__ MX(x),      K__ MX,(x)=K__ MX(x)/   for x,
     K__ MX(x)   for x,
       __ E,(x,α,e)=__ E(x,α,e),   
     __ E(x,α,e),    κ_(x)  = κ(x)        for x,
     κ(x)        for x,
    
       η_(x,α)  = η(x,α)       for x,
    
    η(x,α)       for x,

while the other visco-elastic moduli K__ E, K__ KV, and 
G__ KV are kept fixed. The following statement will be made more specific
and proved in the Appendix:


Let be some solution (u_,_,α_,ϕ_) 
of the initial-boundary-value problem system+–BC-IC 
with the data G-K-eps, which do exist due to Proposition <ref>.
Then (u_,_|_^,α_|_^,ϕ_) 
converge (in terms of subsequences) for →0 to weak solutions to 
system++ with the initial/boundary conditions 
(<ref>a,b)–BC-IC-alpha+.
In particular, a weak solution to the initial-boundary-value problem
system++–(<ref>a,b)–BC-IC-alpha+ 
in the sense of Definition <ref> below does exist.
In addition, this solution conserves energy.

Moreover, the energy dissipated through the Maxwell-type attenuation and 
by damage in
the fluidic regions  over the time interval I converges to zero, i.e. 

    ∫_3/2K__ MX,| sph _|^2
    +G__ MX,| dev _|^2
    +α_η(α_) x̣ṭ→0    for  →0.



It should be emphasized that, although Max-limit similarly as KV-limit
below is intuitively quite clear and generally expected in geophysical 
literature,
its rigorous proof is not entirely trivial and relies on a regularity
which is rather automatic in linear systems but may be highly nontrivial or 
even false in nonlinear hyperbolic systems, as basically also here, 
cf. <cit.>.


Instead of G__ MX→∞ and G__ E(α)→0,
one can think about sending G__ KV→∞. Then, assuming 
the initial condition dev e_ el|_t=0=0, we would have
 dev e_ el→0. Hence, in the limit we would see the model
from Fig. <ref> only with G__ KV replaced
by G__ MX. Since the viscosity in fluidic regions is typically small
(as discuss in particular in the next section <ref>) 
like also G__ KV while G__ MX is typically large in
solid-like geophysical materials, our choice G-K-eps is
more straightforward.
Another difference would be that, instead of π→0,
we would have only sph π→0 while 
dev π→ dev e(u).



Actually, the limit K__ MX→∞ is relevant also in the solid
part . This causes sph π→0 and in the limit one
obtains a model where π is trace-free and accumulates only
shear strain, which is the most common ansatz in plasticity and creep, and in 
geophysical modelling, too.



Passing simultaneously K__ E→∞, we obtain in the limit the Stokes 
incompressible fluid. The bulk viscosity K__ KV then becomes irrelevant, 
so that only the shear viscosity G__ KV remains relevant in 
Figure <ref>. Such model is a great idealization and
has limited application because, beside v__ S→0, now nonphysically 
v__ P→∞. In particular, in this limit, sph e(u) is 
constant in time, namely sph e(u(x,t))=1/3 div u_0(x).
Therefore div  u=0, which expresses the usual incompressibility
condition. Such models are often used in geophysics in short-time-scale
models when the convective term ϱ( u·∇) u, 
occurring in Navier-Stokes equation, can be neglected. 
Cf. e.g. <cit.> for self-gravitating incompressible Stokes model
in layered geoid.




§ SUPPRESSING VISCOSITY TOWARDS PURELY ELASTIC FLUID



As already said, the viscosity (K__ KV,G__ KV) 
in the fluidic domains is only very small (and even not certainly known 
in deep parts of the outer core). Note also that, when G__ E=0
in the fluidic domains, G__ KV is in the position of the 
Maxwellian viscosity. For example, in mantle, 10^2-10^4Pa s is 
considered in <cit.>.
As for the fluidic regions, as said above, the oceans exhibit viscosity around 
10^-3Pa s while the outer core around 10^-2-10^4Pa s.
In any case, these viscosities are much smaller than
typical viscosity in the crust which is of the order 10^22-10^24Pa s
or in the inner core of the order about 10^14-10^15Pa s, 
see e.g. <cit.>.
This certainly gives a good motivation to study the asymptotics when this 
viscosity goes to zero.
In the limit, it yields an inviscid, purely elastic model where the 
Hooke elasticity counts only with the spherical response while the 
shear-stress-free response imitates the ideal inviscid fluid;
sometimes, such sort of models are called compressible Euler fluids.
Such materials are called elastic fluids,
or more specifically just compressible inviscid fluids.
These fluids allow still for propagation of 
P-wavesP-wave!in elastic fluids 
(whose speed is ϱ^-1/2K__ E^1/2 as already mentioned) 
while S-waves are completely excluded.

The resulted system is again (<ref>a,d-f) but now with 

    σ_ sph=3K__ KV sph(e( u)-π)
    +3K__ E sph(e(u)-π)         in ,
    
    3K__ E sph e(u)   in ,
         
    σ_ dev=
    2G__ KV dev(e( u)-π)
    +
    ∂_e^__ E(α, dev(e(u)-π))
       in ,
    
    0   in .


The weak formulation for the limit problem as far as the 
force equilibrium 
(<ref>a,d) 
with system+++ 

concerns 
arises when testing system-u++ by a smooth test function 
v with v(T)=0 and dev e(v)=0 on  and making 
by part integration in time and applying Green formula in space, cf.Definition <ref> in the Appendix.

We can justify this limit-passage scenario rigorously. To this goal, let us 
choose 

    K__ KV,(x)= K__ KV(x),   
     K__ KV(x),      
    G__ KV,(x)= G__ KV(x)   for x,
    
     G__ KV(x)   for x,

while G__ E(α)=0 and K__ MX=G__ MX=∞ on , resulting
already from the limit in Section <ref>.
The following statement
will be proved in the Appendix:


If the initial conditions are enough smooth (cf. the assumption in
Definition <ref>

below),
the solutions (u_,_,α_,ϕ_) 
of the initial-boundary-value problem 
system++–(<ref>a,b)–BC-IC-alpha+
with the data KV-eps, which does exist due to Proposition <ref>,
converge (in terms of subsequences) for →0 to weak solutions to 

(<ref>a,d-f)–(<ref>a,b)–BC-IC-alpha+
with system+++.
In particular, a weak solution to this initial-boundary-value problem


in the sense of Definition <ref> below does exist. Moreover,
this solution conserves energy

and the energy dissipated through the Kelvin-Voigt-type attenuation in
the fluidic regions  over the time interval I converges to zero, i.e. 

    ∫_3/2K__ KV| sph e( u_)|^2
    +G__ KV| dev e( u_)|^2 x̣ṭ→0
        for  →0.







§ APPENDIX: ANALYSIS SKETCHED



In what follows, we will use the (standard) notation for the Lebesgue 
L^p-spaces and W^k,p for Sobolev spaces whose k-th distributional 
derivatives are in L^p-spaces, and the abbreviation H^k=W^k,2. 
Moreover, we will use the standard notation p'=p/(p-1).








In the vectorial case, we will write L^p(Ω;^3)≅ L^p(Ω)^3 
and W^1,p(Ω;^3)≅ W^1,p(Ω)^3. 




For the fixed time interval I=[0,T], we denote by L^p(I;X) the 
standard Bochner space of Bochner-measurable mappings I→ X with X a 
Banach space. Also, W^k,p(I;X) denotes the Banach space of mappings 
from L^p(I;X) whose k-th distributional derivative in time is also in 
L^p(I;X). The dual space to X will be denoted by X^*.



The scalar product between vectors, matrices, or 3rd-order tensors
will be denoted by “·”, “:”, or “”, 
respectively. Finally, in what follows, K denotes a positive, possibly 
large constant.



We will impose the basic assumptions on the data and incorporate them
directly into Definitions <ref>-<ref>. In particular,

    inf_x∈Ω^ϱ(x)>0,    ϱ=0  on ^3∖Ω,    ϱ∈ L^∞(Ω)∩ W^1,3(Ω),
    
       ϱ_ ext W^1,1(I;L^∞(^3)),  ϱ_ ext=0
     outside a bounded set in I×^3,
    
    
       
    K__ E,G__ KV,K__ KV,G__ MX,K__ MX,__ E(·,α,e),
    κ:Ω→[0,∞)  measurable and
               inf_x∈Ωmin(G__ KV(x),K__ KV(x),G__ MX(x),K__ MX(x),K__ E(x),κ)
    >0,   
               sup_x∈max(K__ KV(x)/K__ MX(x),G__ KV(x)/G__ MX(x))<1,
    
    
       __ E(x,·,·):
    [0,1]×^3×3→[0,∞)  twice continuously differentiable, 
               __ E(x,α,·):^3×3→[0,∞)  convex,
    
               
    |_e__ E(x,α,e)|≤ C(1+|e|),   |_α__ E(x,α,e)|≤ C(1+|e|),   
               
    |_e e^2__ E(x,α,e)|≤ C,     |_α e^2__ E(x,α,e)|≤ C,
    
       η:→   uniformly convex, η(0)=0,
    
       u_0 H^1(Ω;^3),  
    v_0 L^2(Ω;^3),  π_0 L^2(Ω;^3× 3),  π_0=0  on  ,
    
                    and α_0 H^1(Ω)  with 0≤α_0(x)≤1 for all xΩ.




Let us emphasize that ass-small-KV/MX is
used only for the second limit passage for estimation 
est-regular–est-regular+ and
is well satisfied in geophysical models where the ratio
of the Kelvin-Voigt and Maxwell viscosities in solid regions is surely
below 10^-8, 
as said in Section <ref>. Also, let us emphasize that
the growth restrictions imposed on 
_α^__ E and 
_α e^2__ E in ass-growth-G
are compatible with the ansatz KV-limit+ and used 
in system-alpha-w and for est-regular below.


We further

integrate system-alpha-w- over Q and apply the Green theorem and the by-part
integration in time to the term (v-α) div(κ∇α),
cf. Remark <ref>. 













The quadruple (u,π,α,ϕ) is called a weak solution to 
the initial-boundary-value problem system+–BC-IC 
provided the data satisfies ass and 

    u W^1,∞(I;L^2(Ω;^3))∩ H^1(I;H^1(Ω;^3)),
    
       π H^1(I;L^2(Ω;^3× 3)),
    
       α∈ L^∞(I;H^1(Ω))∩ H^1(I;L^2(Ω)),
    
       ϕ L^∞(I;H^1(^3)),
 
and the integral identity 

    ∫_Qσ:e(v)-ϱ u· v+f·v x̣ṭ=∫_Ωϱ v_0·v(0) x̣holds for any v∈ H^1(Q;^3) with v(T)=0 and with σ and f from system-u+,
    
    
    
    
    
    
    
    
    
    
       ∫_Q(v-α)
    ∂_α__ E(α, dev(e(u)-π))
    
    +κ∇α∇ v+η(v) x̣ṭ
                   
    +∫_Ωκ/2|∇α_0|^2 x̣≥∫_Qη(α) x̣ṭ
    +∫_Ωκ/2|∇α(T)|^2 x̣for all v∈ L^2(I;H^1(Ω)), 
    
       ∫_I×^3(1/4π g∇ϕ+ϱ u)·∇ v+
    (ϱ+ϱ_ ext)v x̣ṭ=0
 
for all v∈ L^2(I;H^1(^3)), 

and also system-pi holds a.e. on Q, and eventually the resting 
initial conditions u(0)=u_0, π(0)=π_0, and α(0)=α_0 hold 
a.e. on Ω.



Sketched proof of Proposition <ref>.
Without going into details, we may expect that we applied some approximation
method (e.g. a Galerkin-type approximation in space) to obtain some
approximate solution to the 
initial-boundary-value problem system+–BC-IC
which, in principle, can be implemented on 
computers e.g. by the finite-element method. Then the a-priori 
estimates in the space as in est hold also for the approximate
solutions, which may be interpreted as a numerical stability 
of the specific approximation scheme. This approximation leads (after 
smoothening of the potential η) to the initial-value problem for a
system of ordinary-differential equations. Existence of its solution, let
us denote it by (u_k,π_k,α_k,ϕ_k) with k∈ referring to the 
finite-dimensional subspaces used for the Galerkin discretisation, can be 
proved by successive-continuation argument, using the L^∞(I)-estimates 
below. Here we may assume that the initial conditions lie in the 
finite-dimensional spaces used for the Galerkin approximation so that no
further approximation is needed. 

The energy balance engr can serve to see basic apriori estimates and
to perform analysis of the system system+.
Integrating engr over [0,t] and using the by-part integration in time 
for the power of tidal load ϱ_ extϕ and the notation 
(<ref>b,c), we write

    ∫_Ω(ϱ/2| u_k(t)|^2+
    
    3/2K__ E| sph
     e_ el,k(t)|^2
    +__ E(α_k(t), dev e_ el,k(t))
    
    
    
    
       
    +κ/2|∇α_k(t)|^2
    +
    
    1/2ϱ|ω·(x+u_k(t))|^2) x̣
    +∫_^31/
    8 g|∇ϕ_k(t)|^2 x̣
    
       
    +∫_0^t∫_Ω__ KV  e_ el,k∷ e_ el,k
    +__ MXπ_k∷π_k+α_kη'(α_k) x̣ṭ
    
       ≤∫_Ω
    -ϱ∇ϕ_k(t) u_k(t)-ϱϕ_k(t) x̣
       +∫_^3ϱ_ ext(t)ϕ_k(t) x̣
    -∫_0^t∫_^3ϱ_ extϕ_k x̣ṭ+E_0

with e_ el,k=e(u_k)-π_k and with 
the upper bound for the initial energy 

    E_0   =∫_Ωϱ/2|v_0|^2+
    
    3/2K__ E| sph
     e_ el,0|^2
    +__ E(α_0,e_ el,0)
    
    
    
    +ϱϕ(0)+ϱ∇ϕ(0) u_0
    
              +1/2ϱ|ω·(x+u_0)|^2
    
    +κ/2|∇α_0|^2 x̣
    +∫_^3ϱ_ ext(0,x)ϕ(0,x) x̣,

where e_ el,0=e(u_0)-π_0 and where
ϕ(0,x)∈ H^1(Ω) is the gravitational potential 
solving system4 for t=0. The inequality in engr-int have arisen 
from the energy equality engr written for the approximate solution 
and integrated over the time interval [0,t] by forgetting 
some parts of the centrifugal potential with a guaranteed sign.

The non-coercive contribution of the centrifugal force on
the right-hand side of engr-int is to be estimated by Hölder's inequality
as 

    ∫_Ω1/2ϱ|ω|^2|x+u_k(t)|^2
    
     x̣   =∫_Ω1/2ϱ|ω|^2
    |x+u_0(x)+∫_0^t u_k(τ,x) τ̣|^2 x̣
       ≤ C(1+∫_0^t u_k(τ)_L^2(Ω;^3)^2 τ̣)

with some constant C, and then treated by Gronwall's inequality, exploiting 
the kinetic energy on the left-hand side of engr-int.
Furthermore, we can estimate ∫_Ω∇ϕ_k·u_k x̣≤ϵ∇ϕ_k_L^2(Ω;^3)^2
+C_ϵ(∫_0^t u_k_L^2(Ω;^3)^2 ṭ+u_0_L^2(Ω;^3)^2).
Also use the estimate ϕ_k_L^2(Ω)≤ C∇ϕ_k_L^2(^3;^3) 
relying on the boundedness of Ω, provided ϕ(∞)=0  
which is a standard “boundary”
condition for the gravitational potential used in geophysics.





The last integral on the right-hand side of engr-int bears the estimation



    ∫_0^t∫_^3ϱ_ extϕ_k x̣ṭ≤∫_0^tϱ_ ext_L^∞(^3)ϕ_k_L^1(^3)ṭ
       ≤ C∫_0^tϱ_ ext_L^∞(^3)∇ϕ_k_L^2(^3;^3)ṭ≤∫_0^tϱ_ ext_L^∞(^3)(1+ϕ_k_H^1(^3)^2) ṭ,

where we used ϕ_k_
L^1(^3)≤ C∇ϕ_k_L^2(^3;^3);

here the “boundary” condition ϕ_k(∞)=0 together with 
ϱ+ϱ_ ext compactly supported is used.
Altogether, using the Gronwall inequality, we obtain the a-priori estimates 
for the approximate solution:

    u_k_W^1,∞(I;L^2(Ω;^3))∩ 
    H^1(I;H^1(Ω;^3))^≤ C,
    
       π_k_H^1(I;L^2(Ω;^3× 3))^≤ C,
    
       α_k_L^∞(I;H^1(Ω))∩ H^1(I;L^2(Ω))^≤ C,
    
       ϕ_k_L^∞(I;H^1(^3))^≤ C.
 
By Banach's selection principle, we consider a weakly* convergent subsequence
respecting the topologies specified in est+. 
For the limit passage in the nonlinear term _α__ E in 
system-alpha-w, we need to improve it for the 
strong convergence of e_ el,k→ e_ el in L^2(Q;^3×3). 
To this goal, we use the test function u_k-u and π_k-π for 
the Galerkin approximation of system-u+ and system-pi, respectively,
integrate it over the time interval [0,t], and estimate

    ∫_Ω1/2__ KV(e_ el,k(t)-e_ el(t))∷(e_ el,k(t)-e_ el(t))
    
    x̣
       ≤1/2∫_Ω__ KV(e_ el,k(t)-e_ el(t))∷(e_ el,k(t)-e_ el(t))
    +__ MX(π_k(t)-π(t))∷(π_k(t)-π(t))
    
    x̣
    
    
    
    
    
    
       
    +∫_0^t∫_ΩK__ E| sph(e_ el,k-e_ el)|^2
    +(_e^__ E(α_k,e_ el,k)
    -_e^__ E(α_k,e_ el))
    ∷(e_ el,k-e_ el) x̣ṭ
    
    
    
    
    
       =∫_0^t∫_Ω-(ϱ u_k+f_k)(u_k-u)
    
    
           -(__ KV e_ el+K__ E sphe_ el
    +_e^__ E(α_k,e_ el))∷(e_ el,k-e_ el)
    -__ MXπ∷(π_k-π) x̣ṭ
       =
    
    ∫_0^t∫_Ωϱ u_k( u_k- u)-f_k(u_k-u)
    +(_e^__ E(α,e_ el)-_e^__ E(α_k,e_ el)
    )∷(e_ el,k-e_ el)
    
            -(__ KV e_ el+K__ E sphe_ el
    +_e^__ E(α,e_ el))∷(e_ el,k-e_ el)
    -__ MXπ∷(π_k-π) x̣ṭ
                                       
    -∫_Ωϱ u_k(t)(u_k(t)-u(t)) x̣ → 0

with f_k=ϱ(ω×(ω×(x+u_k))+
2ω× u_k+∇ϕ_k); note that f_k is 
bounded in L^2(Q;^3) due to the a-priori estimates (<ref>a,d).
Actually, KV-damage-strong-e(t) is to be understood rather as 
a conceptual strategy:
the mentioned test functions u_k-u and π_k-π
are not legitimate for the Galerkin approximation and (u,π) is still
to be approximated strongly to be valued in the respective finite-dimensional
subspaces - we omitted these standard technical details for simplicity. 
For the convergence to 0 in KV-damage-strong-e(t), we used that 
u_k→ u strongly L^2(Q;^3) due to
the Aubin-Lions theorem, relying on the estimate est-u+ together
with an information about u from the equation system-u+ itself,
and also u_k(t) is bounded in L^2(Ω;^3)
while  u_k(t)→ u(t) strongly in L^2(Ω;^3)
due to the Rellich compact embedding H^1(Ω)⊂ L^2(Ω),




and also 
_e^__ E(α_k,e_ el)→_e^__ E(α,e_ el)
strongly in L^2(Q;) 
due to the continuity of the Nemytskiĭ mapping 
induced by _e^__ E(·,e_ el) and 
α_k→α in L^2(Q) again just by the Rellich theorem
since both α_k and ∇α_k is estimated in L^2-spaces.

Thus, from KV-damage-strong-e(t), we obtain
e_ el,k(t)→ e_ el(t) strongly in L^2(Ω;^3× 3) 
for all t∈ I. Using 
it for a general t I, we obtain e_ el,k→ e_ el strongly in 
L^2(Q;^3× 3) by the Lebesgue theorem.








The limit passage in the Galerkin approximation towards the 
integral identities system-w is then simple by weak/strong
continuity or, in case of the variational inequality system-alpha-w, 
also semicontinuity.

For the energy conservation, the essential needed facts are that 
√(ϱ) u∈ L^2(I;H^1(Ω;^3)^*) is in duality with 
√(ϱ) u∈ L^2(I;H^1(Ω;^3)) and also that 
div(κ∇α)∈ L^2(Q) is in duality with 
α∈ L^2(Q) so that the by-part
integration formulas rigorously hold:

    ∫_0^t⟨√(ϱ) u,√(ϱ) u⟩ ṭ
    =∫_Ωϱ/2| u(t)|^2-ϱ/2| u(0)|^2 x̣,    and
       ∫_0^t∫_Ω div(κ∇α)α x̣ṭ
    =∫_Ωκ/2|∇α(0)|^2
    -κ/2|∇α(t)|^2 x̣ ;
 
see e.g. <cit.>. More in detail, for 
calculus1 we have used the comparison 
√(ϱ) u=( div σ-f)/√(ϱ)

and the estimate 

    √(ϱ) u_L^2(I;H^1(Ω;^3)^*)
    =
    sup_v_L^2(I;H^1(Ω;^3))≤1∫_Q div σ-f/√(ϱ)v x̣ṭ
               
    =sup_v_L^2(I;H^1(Ω;^3))≤1∫_Qσ∷(√(ϱ)∇ v-
    v⊗∇ϱ/√(ϱ))
    +f/√(ϱ)v x̣ṭ≤ C

for which the smoothness ϱ∈ W^1,3(Ω) is needed due to
the occurrence of ∇ϱ in est-of-DDTu, 
while for calculus2
we have used the comparison div(κ∇α)∈ζ(α)
+_α__ E(α, dev e_ el).




The quadruple (u,π,α,ϕ) is called a weak solution to 
the initial-boundary-value problem  
(<ref>a,b)–system++–BC-IC-alpha+
with the data G-K-eps provided again the data satisfies ass and 
(<ref>a,d) holds together with 

    π H^1(I;L^2(;))
       and   α∈ L^∞(I;H^1())∩ H^1(I;L^2())

and the integral identity momentum-w holds with 
σ=σ_ sph+σ_ dev from (<ref>b,c), 
also system-alpha-w holds with 
Q and Ω replaced respectively  and ,
furthermore gravity-w holds, 
and also system-pi holds a.e. on , and eventually 
u(0)=u_0 holds a.e. on Ω, and π(0)=π_0 and 
α(0)=α_0 hold a.e. on .



Sketched proof of Proposition <ref>.

Like engr-int but now with G-K-eps taken into account, we have

    ∫_Ωϱ/2| u_(t)|^2
    +3/2K__ E| sph e_ el,(t)|^2+
    1/2ϱ|ω·(x+u_(t))|^2 x̣
    +∫_^31/
    8 g|∇ϕ_(t)|^2 x̣
       
    +∫___ E(α_(t), dev e_ el,(t))
    +κ/2|∇α_(t)|^2 x̣
       +
    ∫___ E(α_(t), dev e_ el,(t))
    +κ/2|∇α_(t)|^2 x̣
    
    
    
       
    +∫_0^t∫_Ω__ KV  e_ el,∷ e_ el, x̣ṭ
    +∫_0^t∫___ MXπ_∷π_
    +α_η'(α_) x̣ṭ
       
    +∫_0^t∫_1/__ MXπ_∷π_
    +α_η'(α_) x̣ṭ
    ≤∫_Ω
    -ϱ∇ϕ_(t) u_(t)-ϱϕ_(t) x̣
       +∫_^3ϱ_ ext(t)ϕ_(t) x̣
    -∫_0^t∫_^3ϱ_ extϕ_ x̣ṭ+E_0

with E_0 again from def-E0, relying that ≤1.
By the Gronwall-inequality arguments like used for engr-int,
we obtain the a-priori estimates (<ref>b,d) now for (π_,ϕ_)
instead of (π_k,ϕ_k) and also 

    u__W^1,∞(I;L^2(Ω;^3))∩ 
    H^1(I;H^1(;^3))^≤ C
    
       α__L^∞(I;H^1())∩ H^1(I;L^2())^≤ C,
    
       π__L^2(;^3× 3)^≤√()C,
    
       α__L^2(;^3× 3)^≤ C/√(),
        and    
       ∇α__L^∞(I;L^2(;^3))≤ C/√().

Therefore, for a subsequence, we have

    u_→ u      weakly* in W^1,∞(I;L^2(Ω;^3))∩ 
    H^1(I;H^1(Ω;^3)),
    
       π_→π      weakly in H^1(I;L^2(Ω;)),
    
       π_|_^→0      strongly in  H^1(I;L^2(;)),
    
       α_|_^→α      weakly* in L^∞(I;H^1())∩ H^1(I;L^2()),
    
       ϕ_→ϕ      weakly* in  L^∞(I;H^1(^3)),  and also
       
    e_ el,→ e_ el      strongly in L^2(Q;).
 
For conv-pi-L, we used π_→0 strongly in 
L^2(;^3× 3) due to est-pi-2 together with the assumption 
π_0=0 on .
For conv-el-2, like KV-damage-strong-e(t) when taking into
account that π=0 on , we have

    ∫_Ω1/2__ KV(e_ el,(t)-e_ el(t))∷(e_ el,(t)-e_ el(t)) 
    x̣
       ≤∫_Ω1/2__ KV(e_ el,(t)-e_ el(t))∷(e_ el,(t)-e_ el(t)) x̣
               
    +∫_1/2__ MX(π_(t)-π(t))∷(π_(t)-π(t))
    
     x̣+∫_1/2__ MXπ_(t)∷π_(t) x̣
    
    
    
    
    
    
       
    +∫_0^t∫_ΩK__ E| sph(e_ el,-e_ el)|^2
    +(_e^__ E,(α_,e_ el,)
    -_e^__ E,(α_,e_ el))
    ∷(e_ el,-e_ el) x̣ṭ
    
    
    
    
    
       =∫_0^t(∫_Ω-(ϱ u_+f_)(u_-u)
    
    -(__ KV e_ el+K__ E sphe_ el
    +_e^__ E,(α_,e_ el))∷
                                   ∷(e_ el,-e_ el) x̣
    -∫___ MXπ∷(π_-π) x̣) ṭ
       =
    
    ∫_0^t(∫_Ωϱ u_( u_- u)
    -f_(u_-u)
    +(_e^__ E,(α,e_ el)-_e^__ E,(α_,e_ el)
    
    
                     -
    
    __ KV e_ el-K__ E sphe_ el
    -_e^__ E,(α,e_ el))∷(e_ el,-e_ el) x̣
               
    -∫___ MXπ∷(π_-π) x̣) ṭ
    -∫_Ωϱ u_(t)(u_(t)-u(t)) x̣ → 0

with f_=ϱ(ω×(ω×(x+u_))+
2ω× u_+∇ϕ_)
bounded in L^2(Q;^3) and where we now used 
the continuity of the Nemytskiĭ mapping induced by 
_e^__ E,(·,e_ el)=_e^__ E(·,e_ el) 
on  while, on , we use that simply 
_e^__ E,(α_,e_ el)=
_e^__ E(α_,e_ el)→0 strongly in 
L^2(); recall the scaling G-K-eps.


The limit passage towards the weak solution in the sense of 
Definition <ref> is then simple by the continuity with 
respect to the convergences conv-2. More details deserve 
only the limit passage in 
system-alpha-w which, written for (u_,π_,α_) and
omitting the terms ∫_η(α_) x̣ṭ≥0 and
∫_κ/2|∇α_(T)|^2 x̣≥0,
reads as:

    ∫_(v-α_)
    ∂_α__ E,(α_, dev(e(u_)-π_))
    
    +κ∇α_∇ v+η(v) x̣ṭ
    +∫_κ/2|∇α_0|^2 x̣
           
    +(∫_(v-α_)
    ∂_α__ E,(α_, dev(e(u_)-π_))
    +κ∇α_∇ v+η(v) x̣ṭ
                 
    +∫_κ/2|∇α_0|^2 x̣)
    ≥∫_η(α_) x̣ṭ
    +∫_κ/2|∇α_(T)|^2 x̣.



Now we use 

    |∫_α_∂_α__ E,(α_,
    
    e_ el,)
    x̣ṭ|
    ≤ Cα__L^2()^1+|e_ el,| _L^2(Q)^
    =𝒪(√())→0
 
due to est-DT-a-2, and

    |∫_κ∇α_∇ vx̣ṭ|
    ≤(sup_Ωκ))∇α__L^2(;^3)^∇ v_L^2(Q;^3)^
    =𝒪(√())→0
 
due to est-nabla-a-2.
In the limit we thus obtain system-alpha-w on .

Altogether, we proved that the limit (u,π,α,ϕ)
solves the initial-boundary-value problem
system++–(<ref>a,b)–BC-IC-alpha+
in the sense of Definition <ref>.

In addition, this solution conserves energy. This can be shown
again by using that √(ϱ) u in duality with 
√(ϱ) u so that calculus1 holds, and 
that div(κ∇α)∈ L^2() and 
α∈ L^2() so that also calculus2 holds
but now with  and  instead of Q and Ω, respectively.

To prove Max-limit, we use engr-int written for 
(u_,π_,α_,ϕ_) with t=T and the scaling 
BC-IC-alpha+, namely

    lim sup_→0∫_1/__ MXπ_∷π_
    +α_ζ(α_) x̣ṭ
       ≤lim_→0∫_^3ϱ_ ext(T)ϕ_(T) x̣
    -∫_Ωϱ∇ϕ_(T) u_(T)-ϱϕ_(T)
    
     x̣
    -∫_0^T∫_^3ϱ_ extϕ_ x̣ṭ
    
       
    -∫_Q(ω×(ω×(x+u_)· u_ x̣ṭ
    -∫___ E(α_(T),e_ el,(T)) x̣
       
    -lim inf_→0(∫_Ωϱ/2| u_(T)|^2
    
    +3/2 K__ E| sph e_ el,(T)|^2
    
    
     x̣
       +∫___ E(α_(T),e_ el,(T))
    +κ/2|∇α_(T)|^2 x̣
    +
    ∫_^31/
    8 g|∇ϕ_(T)|^2 x̣
       
    +∫_Q__ KV e_ el,∷ e_ el,
     x̣ṭ+∫___ MXπ_∷π_
    +α_η'(α_) x̣ṭ)
    +E_0
    
       ≤∫_^3ϱ_ ext(T)ϕ(T)-1/
    8 g|∇ϕ(T)|^2 x̣
    -∫_Ωϱ∇ϕ(T) u(T)-ϱϕ(T)
    
     x̣
       
    -∫_0^T∫_^3ϱ_ extϕ x̣ṭ
    -∫_Q__ KV e_ el∷ e_ el+(ω×(ω×(x+u)· u x̣ṭ
    
       
    -∫_Ωϱ/2| u(T)|^2
    
    +3/2 K__ E| sph e_ el(T)|^2
    
    
     x̣
    -∫___ E(α(T),e_ el(T))+κ/2|∇α(T)|^2 x̣
    
       
    
    
    
    
    -∫___ MXπ∷π
    +αη'(α) x̣ṭ+E_0=0,

where we again used the notation C-D and E_0 from 
def-E0 now with κ=0 on . The 
first inequality arises from “forgetting” the 
term κ/2|∇α_(T)|^2 on , while the
second inequality is by weak lower semicontinuity.
The last equality in engr-int-eps expresses the energy conservation 
for the limit system, discussed already above. 





















Thus Max-limit is proved.


From KV-damage-strong-e(t)-2, one can even see that also π_→π
strongly so that, together with conv-el-2, even the total strain
e(u_)=e_ el,+π_ converges strongly in L^2(Q;). We 
did not need this additional result in the above convergence proof, however.



Since ∂ζ(α) is bounded in L^2(Q) in our model,
by comparison, we have also 
div(κ∇α)∈η(α)
+∂_α__ E(α, dev e_ el) bounded in L^2(Q), cf.system-alpha++. Since also α L^2(Q), 
the formula ∫_Qα div(κ∇α) x̣ṭ
=1/2∫_Ωκ|∇α_0|^2-κ|∇α(T)|^2 x̣
rigorously holds, and we can write system-alpha-w in
as the original inequality system-alpha-w- holding a.e. on Q.
The integral form system-alpha-w is however suitable for the limit
passages, in contrast to system-alpha-w-.










Assuming, beside ass, also 
e_ el,0=e(u_0)-π_0∈ H^1(Ω;^3×3) and 
v_0∈ H^2(Ω;^3), then the weak solution to the system 
(<ref>a,d-f)–(<ref>a,b)–BC-IC-alpha+
with system+++ is understood as a five-tuple
(u,σ_ sph,σ_ dev,π,ϕ)∈ W^1,∞(I;L^2(Ω;^3))× L^2(Q;)× L^2(Q;)× H^1(I;L^2(;))× W^1,∞(I;H^1(^3))







if the integral identity

    ∫_σ_ dev∷ dev e(v) x̣ṭ+
    ∫_Q
    σ_ sph∷ sph e(v)-uv+f v x̣ṭ
    =∫_Ωv̊_0 v(0,·) x̣

with f from system-u+ holds for any v H^1(Q;^3) with v(T)=0, 

further system-pi++ hold a.e. on  and system+++ relating 
(σ_ dev,σ_ sph) with (u,π) hold a.e. on Q, and also 
the initial conditions u(0,·)=u_0 and π(0,·)=π_0 hold 
a.e. on Ω.


Note that, controlling σ_ sph and σ_ dev in the 
above definition, we have implicitly also included the information 

    u|_ H^1(I;H^1(Ω;^3))    and     
    u L^2(I;L_ div^2(Ω;^3))

with L_ div^2(Ω;^3)
={u L^2(Ω;^3);  div u L^2(Ω)},

while the strain e(u) is not defined in the fluidic regions .
Also one has e(u)-π H^1(I;L^2(;)).




Sketched proof of Proposition <ref>.
We first prove a certain regularity by differentiating in time the system
(<ref>a,d) with system+++ written for the solution obtained in 
Proposition <ref> and denoted now as (u_,π_,α_,ϕ_), 
and employing the test by ( u_,π_). Taking into account the 
scaling KV-eps and using again the orthogonality for the Coriolis force 
(now in terms of accelerations) as (2ϱω× u)· u=0, 
this gives

    ∫_Ωϱ/2| u_(t)|^2
    +3/2K__ E| sph  e_ el,(t)|^2 x̣
           
    +∫_0^t∫___ KV e_ el,∷ e_ el,
    +__ MXπ_∷π_ x̣ṭ
    +∫_0^t∫___ KVe( u_)∷ e( u_) x̣ṭ
       
    =∫_Ωϱ/2| u_(0)|^2
    +3/2K__ E| sph  e_ el,(0)|^2x̣
    -∫_0^t∫_Ωϱ(ω×(ω×(x+ u_))+∇ϕ_) u_
           
    -∫_0^t∫_(_ee^2__ E(α_,e_ el,)∷ dev  e_ el,+
    _eα^2__ E(α_,e_ el,)α_)
    ∷ dev  e_ el,x̣ṭ
    
       ≤1/2ϱ_L^∞(Ω)^ u_(0)_L^2(Ω;^3)^2
    +3/2K__ E sph(e(v_0)-π_0)_L^2(Ω;^3× 3)^2
    
         
    +√(ϱ)_L^∞(Ω)^ω×(ω×(x+ u_))+∇ϕ__L^2(Q;^3)^2+∫_0^t∫_Ωϱ/2| u_(t)|^2 x̣ṭ
         
    +1/c(
    _ee^2__ E(α_,e_ el,)_L^∞(;^3^4)^2
    e( u_)_L^2(Q;^3× 3)^2
    
       
      +
    _eα^2__ E(α_,e_ el,)_L^∞(Q;^3× 3)^2
    α__L^2(Q)^2)
    
    +ϵ
     e_ el,_L^2([0,t]×;^3× 3)^2

with c:=1/2__ KV_L^∞(Ω;^3^4); 
then the last term in est-regular can be absorbed in the left
hand side. Note that we needed the assumptions ass-growth-G to have
_eα^2__ E(α_,e_ el,) and
_ee^2__ E(α_,e_ el,) apriori bounded.
We further use

    u_(0)=1/ϱ( div σ_0,
    +ϱ(ω×(ω×(x+u_0))+
    2ω×v_0+∇ϕ_0)
    )∈ L^2(Ω;^3)

with the initial stress σ_0,=σ_ sph+σ_ dev 
with σ_ sph and σ_ dev from (<ref>b,c)
with (K__ KV,G__ KV)=(K__ KV,,G__ KV,) from KV-eps,

and with e_ el=e(u_0)-π_0 and e_ el=e(v_0)-π_0,
and ϕ_0∈ H^1(^3) solving 
div(∇ϕ_0/(4 g)+ϱ u_0)=ϱ+ϱ_ ext(0).
Note that indeed u_(0)∈ L^2(Ω;^3)
due to the assumptions 
e_ el,0=e(u_0)-π_0∈ H^1(Ω;^3×3) and 
v_0∈ H^2(Ω;^3) involved in Definition <ref>. 
Also note that π_0:=π(0) is involved in est-regular
and implicitly also in est-regular+, although we do not 
prescribe any initial condition on π. Yet, from system-pi++,
we can read π(0)=__ MX^-1σ_0 on  while 
π(0)=0 on  because we have now __ MX^-1=0 on . 
To estimate σ_0 let us realize that
σ_0_H^1(Ω;^3×3)≤__ KV(e(v_0)-π(0))
+_e_ E(α_0,e_ el,0)_H^1(Ω;^3×3)≤
C+__ KV__ MX^-1σ_0_H^1(Ω;^3×3)
with some C depending on v_0_H^2(Ω;^3), α_0_H^1(Ω),
and e_ el,0_H^1(Ω;^3×3); note that these quantities 
have been supposed finite in Definition <ref>.
Thus, using ass-small-KV/MX, we get the desired bound on σ_0.
Hence, e_ el(0) which occurs in est-regular
itself can be estimated.

Using Gronwall's inequality for est-regular then yields the estimates 

    u__W^2,∞(I;L^2(Ω;^3)) ∩ H^2(I;H^1(;^3))^≤ C,
    
    
       π__L^2(;^3× 3)^≤ C,   and   
       u__H^2(I;H^1(;^3))^≤ C/√(),
    which are now at disposal together with est-phi and 
    est-a-2.
    
    Also we have   σ__L^2(Q;)^≤ C   with  σ_=__ KVe( u_)+3K__ KV sph e(u_)   on ,
      __ KV e_ el,+__ E(α_,e_ el,)
       on ,

where we again use the notation C-D and e_ el,=e(u_)-π_.

By Banach's selection principle, we consider a weakly* convergent subsequence
respecting the topologies specified in (<ref>a,b) and est-a-2
together with 
the convergence conv-phi and also σ_→σ weakly in 
L^2(Q;). Then we put naturally σ_ sph:= sph σ and 
σ_ dev:= dev σ. 
As e( u_)_L^2(;^3×3)=𝒪(1/√()),
we have __ KVe( u_)_L^2(;^3×3)
=𝒪(√())→0 so that
σ=3K__ KV sph e(u) on .
To identify σ=__ KV e_ el+__ E(α_,e_ el)
in the solid regions ,
due to nonlinearities __ E(α,·), we again need to prove 
e_ el,→ e_ el strongly in 
L^2(;). To this goal, one is to modify 
KV-damage-strong-e(t)-2 to be used on  rather than Ω.
The peculiarity is that u_-u is no longer a legitimate
test function because e(u) is not well defined in the fluidic regions 
. For this reason, we take some smooth ũ_ that will
approximate u strongly in L^∞(I;L_ div^2(Ω;^3)) 
and even ũ_|_→ u|_ strongly in 
H^1(I;H^1(Ω;^3)), and we can assume that this convergence is 
sufficiently slow so that 

    ũ__L^2(I;H^1(;^3))≤ 1/√().

Then, denoting ẽ_ el,=e(ũ_)-π,
we have ẽ_ el,(t)→ e_ el(t) for a.a. t I
and we can write

    lim sup_→0∫_1/4__ KV(e_ el,(t)-e_ el(t))∷(e_ el,(t)-e_ el(t)) 
    x̣
       ≤lim sup_→0∫_1/2__ KV(e_ el,(t)-ẽ_ el,(t))∷(ẽ_ el,(t)-e_ el,(t)) x̣
                   
    +lim_→0∫_1/2__ KV(ẽ_ el,(t)-e_ el(t))∷(ẽ_ el,(t)-e_ el(t)) x̣
       ≤lim sup_→0(∫_(1/2__ KV(e_ el,(t)-ẽ_ el,(t))∷(e_ el,(t)-ẽ_ el,(t))
    
    
           
    +
    1/2__ MX(π_(t)-π(t))∷(π_(t)-π(t)))
    
     x̣
    
    
    
    
    
    +∫_0^t∫_ΩK__ E| sph(e_ el,-ẽ_ el,)|^2 x̣ṭ
               
    +∫_0^t∫_(_e^__ E(α_,e_ el,)
    -_e^__ E(α_,ẽ_ el,))
    ∷(e_ el,-ẽ_ el,) x̣ṭ)
    
    
    
    
    
       ≤lim_→0(
    ∫_/2__ KVe(u_0) x̣
    -∫_0^t∫___ KVe( u_)∷ e(ũ_)
     x̣ṭ
    
           
    -
    ∫_0^t∫_Ω(ϱ u_+f_)(u_-ũ_)
    -∫_0^t∫___ MXπ∷(π_-π)
    
       
        
    +(__ KVẽ_ el,
    +K__ E sphẽ_ el,
    +_e^__ E(α_,ẽ_ el,))
    ∷(e_ el,-ẽ_ el,)
     x̣ṭ)=0,

where the third inequality have arisen by “forgetting” the nonnegative
term ∫___ KVe(u_(T))∷ e(u_(T)) x̣.
Note the we used approximation together with est-pi-3+ for the
estimate

    |∫_0^t∫___ KVe( u_)∷ e(ũ_)
     x̣ṭ|≤__ KV_L^∞(Ω;^3^4)×
               ×e( u_)_L^2(;^3×3)e(ũ_)_L^2(;^3×3))
    =𝒪(1/√())
    =𝒪(√())→0.











From KV-damage-strong-e(t)-3, we thus have 
e_ el,|_(t)→ e_ el|_(t) at a.a. t I. Then, 
instead of conv-el-2, by the Lebesgue theorem, we now proved

    e_ el,|_→ e_ el|_      strongly in L^2(;),

which is to be used for the limit passage in the nonlinear term
_e^__ E.

The energy conservation now holds due to the proved regularity, as 
√(ϱ) u∈ L^∞(I;L^2(Ω;^3)) is surely in duality with 
√(ϱ) u∈ W^1,∞(I;L^2(Ω;^3)), cf. est-u-3.

Eventually, like we did in engr-int-eps, we now can show that

    lim sup_→0∫___ KV e( u_)∷ e( u_) x̣ṭ≤∫_^3ϱ_ ext(T)ϕ(T)-1/
    8 g|∇ϕ(T)|^2 x̣
    
       
    -∫_Ωϱ∇ϕ(T) u(T)-ϱϕ(T)
    
     x̣-∫_0^T∫_^3ϱ_ extϕ x̣ṭ
    -∫_Q(ω×(ω×(x+u)· u x̣ṭ
    
       
    -∫_Ωϱ/2| u(T)|^2
    
    +3/2 K__ E| sph e_ el(T)|^2
    
    
     x̣-∫___ E(α(T),e_ el(T))+κ/2|∇α(T)|^2 x̣
    
       
    
    
    
    
    -∫___ KV e_ el∷ e_ el+__ MXπ∷π
    +αη'(α) x̣ṭ+E_0=0,

with E_0 from def-E0. The last equality is due to the mentioned energy 
conservation. Thus KV-limit is proved.



Note that, in est-regular, we benefited from having 
G__ E already pushed to zero on  because the
viscosity on  is (intentionally) not uniformly controlled,
being limited to zero. Thus the direct merging of both 
limit processes in Propositions <ref> 
and <ref> is not possible. Of course, a suitable
scaling between these two would facilitate such connection.


The interfaces between the ocean beds and the mantle as well as the 
Gutenberg's discontinuity and  between the inner and outer core regions 
typically also exhibit discontinuities in
mass density ϱ, which is incompatible with the assumption 
ϱ∈ W^1,3(Ω) in ass-rho used in est-of-DDTu.
Note that the additional regularity of the initial conditions involved
in Definition <ref> allowed us to avoid this restriction and
consider a general, possibly discontinuous ϱ∈ L^∞(Ω).
A respective modification of the proofs of 
Propositions <ref>–<ref> would be possible, too.




Acknowledgments. 
The author is thankful to Katharina Brazda and Ctirad Matyska for
many fruitful discussions about the models. 

Moreover, deep thanks are also for hospitality and support of the 
Erwin Schrödinger Institute, Univ. Vienna, 
and for the partial support of the 
Czech Science Foundation projects 16-03823S and 17-04301S and also 
by the Austrian-Czech projects 16-34894L (FWF/CSF),
as well as through the institutional support RVO: 61388998 (ČR).

99

BeZi01DRRM
Y. Ben-Zion.
Dynamic ruptures in recent models of earthquake faults.
J. Mech. Phys. Solids, 49:2209–2244, 2001.

BZamp09SRRS
Y. Ben-Zion and J.-P. Ampuero.
Seismic radiation from regions sustaining material damage.
Geophys. J. Int., 178:1351–1356, 2009.

Bog77HECV
D. Boger.
A highly elastic constant-viscosity fluid.
J. Non-Newtonian Fluid Mechanics, 3:87–91, 1977.

Braz17EGEG
K. Brazda.
The elastic-gravitational equations in global seismology with
  low regularity.
PhD thesis, Univ. Wien, 2017.

BrHoHo17VFEE
K. Brazda, M. V. de Hoop, and G. Hoermann.
Variational formulation of the earth's elastic-gravitational
  deformations under low regularity conditions.
Preprint arXiv:1702.04741, 2017.

DahTro98TGS
F. A. Dahlen and J. Tromp.
Theoretical global seismology.
Princetown Univ. Press, Princetown, NJ, 1998.

GreNag65GTEP
A. Green and P. Naghdi.
A general theory of an elastic-plastic continuum.
Arch. Rational Mech. Anal., 18:251–281, 1965.

HBAD09DERC
R. A. Harris et al.
The SCEC/USGS dynamic earthquake rupture code verification
  exercise.
Seismological Res. Lett., 80:119–126, 2009.

HuAmHel14ERMW
Y. Huang, J.-P. Ampuero, and D. V. Helmberger.
Earthquake ruptures modulated by waves in damaged fault zones.
J. of Geophysical Research: Solid Earth, B9:3133–3154, 2014.

KaLaAm08SEMS
Y. Kaneko, N. Lapusta, and J.-P. Ampuero.
Spectral element modeling of spontaneous earthquake rupture on rate
  and state faults: Effect of velocity-strengthening friction at shallow
  depths.
J. Geophysical Res., 113:B09317, 2008.

KomTro02SESG
D. Komatitsch and J. Tromp.
Spectral-element simulations of global seismic wave propagation -
  I. validation.
Geophys. J. Int., 149:390–412, 2002.

KomTro02SESG2
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.

KooDum11VEIC
L. Koot and M. Dumberry.
Viscosity of the Earth's inner core: Constraints from nutation
  observations.
Earth and Planetary Science Letters, 308(3):343–349, 2011.

LayWal95MGS
T. Lay and T. C. Wallace.
Modern global seismology.
Acad. Press, San Diego, 1995.

LyaBZ14DBRM
V. Lyakhovsky and Y. Ben-Zion.
Damage-breakage rheology model and solid-granular transition near
  brittle instability.
J. Mech. Phys. Solids, 64:184–197, 2014.

LHAB09NDRW
V. Lyakhovsky, Y. Hamiel, J.-P. Ampuero, and Y. Ben-Zion.
Non-linear damage rheology and wave resonance in rocks.
Geophys. J. Int., 178:910–920, 2009.

LyHaBZ11NLVE
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.

LyaMya84BECS
V. Lyakhovsky and V. P. Myasnikov.
On the behavior of elastic cracked solid.
Phys. Solid Earth, 10:71–75, 1984.

MeaFur13SSWO
T. Maedae and T. Furumura.
FDM simulation of seismic waves, ocean acoustic waves, and tsunamis
  based on tsunami-coupled equations of motion.
Pure Appl. Geophys., 170:109–127, 2013.

PPAB12TDDR
C. Pelties, J. de la Puente, J.-P. Ampuero, G. B. Brietzke, and M. Käser.
Three-dimensional dynamic rupture simulation with a high-order
  discontinuous galerkin method on unstructured tetrahedral meshes.
J. Geophys. Res., 117:B02309, 2012.

RajRou03EDSM
K. R. Rajagopal and T. Roubíček.
On the effect of dissipation in shape-memory alloys.
Nonlinear Anal., Real World Appl., 4:581–597, 2003.

Roub13NPDE
T. Roubíček.
Nonlinear Partial Differential Equations with Applications.
Birkhäuser, Basel, 2nd edition, 2013.

Roub14NRSF
T. Roubíček.
A note about the rate-and-state-dependent friction model in a
  thermodynamical framework of the Biot-type equation.
Geophysical J. Intl., 199:286–295, 2014.

Roub17GMHF
T. Roubíček.
Geophysical models of heat and fluid flow in damageable poro-elastic
  continua.
Cont. Mech. Thermodyn., 29:625–646, 2017.

RoPaMa13QACV
T. Roubíček, C. G. Panagiotopoulos, and V. Mantič.
Quasistatic adhesive contact of visco-elastic bodies and its
  numerical treatment for very small viscosity.
Zeitschrift angew. Math. Mech., 93:823–840, 2013.

RoSoVo13MRLF
T. Roubíček, O. Souček, and R. Vodička.
A model of rupturing lithospheric faults with re-occurring
  earthquakes.
SIAM J. Appl. Math., 73:1460–1488, 2013.

Secc13VOC
R. A. Secco.
Viscosity of the outer core.
In T. Ahrens, editor, Mineral Physics & Crystallography: A
  Handbook of Physical Constants, pages 218–226. Willey, 2013.

SmyPal07VEOC
D. E. Smylie and A. Palmer.
Viscosity of Earth's outer core.
Preprint arXiv:0709.3333, 2007.

ToCaMa09SSLV
N. Tosi, O. Čadek, and Z. Martinec.
Subducted slabs and lateral viscosity variations: effects on the
  long-wavelength geoid.
Geophys. J. Int., 179:813–826, 2009.

TAKS13EEEE
V. Tsai, J.-P. Ampuero, H. Kanamori, and D. Stevenson.
Estimating the effect of Earth elasticity and variable water
  density on tsunami speeds.
Geophysical Res. Letters, 40:492–496, 2013.

DKVD98VLIP
G. A. D. Wijs, G. Kresse, L. Vočadlo, D. Dobson, D. Alfé, M. J. Gillan,
  and G. D. Price.
The viscosity of liquid iron at the physical conditions of the
  Earth's core.
Nature, 392 (6678):805–807, 1998.

WooDeu09TOEF
J. H. Woodhouse and A. Deuss.
Theory and observations - Earth's free oscillations.
In B. Romanowicz and A. Dziewonski, editors, Seismology and
  Structure of the Earth: Treatise on Geophysics, volume 1, chapter 1.02,
  pages 31–65. Elsevier, 2009.



