# Do multiphysics processes lead to mesh independent analyses?

Eleni Gerolymatou, Alexandros Stathas, Ioannis Stefanou

► To cite this version:

Eleni Gerolymatou, Alexandros Stathas, Ioannis Stefanou. Do multiphysics processes lead to mesh independent analyses?. International Journal of Mechanical Sciences, 2024, 274, pp.109265. 10.1016/j.ijmecsci.2024.109265. hal-04934485

HAL Id: hal-04934485

https://nantes-universite.hal.science/hal-04934485v1

Submitted on 21 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Do multiphysics processes lead to mesh independent analyses?

Eleni Gerolymatou<sup>a</sup>, Alexandros Stathas<sup>b,*</sup>, Ioannis Stefanou<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Bifurcation

Solid mechanics

Finite elements

Strain localization

Inelasticity

Multiphysics

## ABSTRACT

Multiphysics processes are common in nature and can occur in conjunction with intense strains. When strain softening is also present, the original uniform strain distribution is lost and strain localization occurs.

In a classical continuum strain localization happens on a mathematical plane (i.e shear band of zero thickness). This leads to mesh dependent numerical analyses, in which strain localizes on the smallest mesh length.

In this paper we investigate under which conditions multiphysics processes regularize the problem of strain localization. We study, for the first time, the influence of all the different coupling terms present in the original coupled multiphysics balance equations on the strain localization profile. We also consider the influence of viscosity and inertia effects.

We show that the multiphysics processes do not always lead to mesh independent results. Furthermore, our rigorous mathematical analysis provides bounds for the different coupling parameters, for which multiphysics processes lead to strain regularization and mesh independent numerical results. Finally, we corroborate our mathematical findings, applying them to a numerical example of a coupled hydro-mechanical system of a fault under shear.

## 1. Introduction

Multiphysics processes are ubiquitous in nature and a large number of them has a significant influence on the mechanical response of materials [1–17]. Common examples are thermal diffusion, water flow and chemical reactions taking place in soil and other composite materials. These processes are of utmost importance during material failure and strain softening.

Considering a material under constant strain, when strain softening is present in the material response, the deformation inside the material may deviate from the homogeneous state leading to strain localization (i.e. strain is no longer homogeneous inside the whole material but localizes instead in smaller regions). From a theoretical point of view, when strain localization occurs, the homogeneous solution characterizing the deformation of the underlying continuum mechanics problem is unstable (see [18,19]). This leads to a zone of intense deformation, whose thickness is related to energy dissipation and post peak strength considerations. When simple continua are taken into account, see [20,21], this zone has zero thickness (localization on a mathematical plane). This is problematic in numerical simulations, because of the analysis results becoming mesh dependent (i.e. the profile of strain localizes to the smallest length in the mesh, without converging as we increase the number of elements, see also [11,19,22–26]).

In order to avoid mesh dependence, different approaches of regularization of strain are in use. These include the use of non-simple continua (see [27–29]), for example the use of micropolar materials in [30] and nonlocal materials in [28]. In these cases the width of the localization zone is controlled by the material properties (i.e. internal length, see [31]). The solution remains non unique in the sense that different strain localization (i.e. perturbation) patterns satisfying the boundary conditions can be developed, but localization now can take place only in zones of finite thickness, avoiding mesh dependency and restoring the dissipative properties of the medium (i.e. strain regularization takes place).

Viscosity in the material response has also been used to a significant extent to resolve the issue of mesh dependence. In this case the width of the localization zone has been shown to be controlled by defects in the initial distribution of material properties (see [32,33]). Thus, the strain localization profile becomes a choice of the analyst and is not connected to an internal length scale characteristic of the problem.

The addition of the inertia was proposed to offer a minimum length-scale for the regularization of the strain localization pattern (see [34–36]). However, it was recently shown in [37], that even in the presence of inertia effects, intrinsic mechanical strain softening, with viscous strain rate hardening does not regularize strain localization on a

<footnote><sup>a</sup>Technische Universität Clausthal, Clausthal-Zellerfeld, Germany</footnote>

<footnote><sup>b</sup>Institut de Recherche en Génie Civil et Mécanique (UMR CNRS 6183), Ecole Centrale de Nantes, Nantes, France</footnote>

<footnote>* Corresponding author.
E-mail address: alexandros.stathas@ec-nantes.fr (A. Stathas).</footnote>
mathematical plane and it does not alleviate mesh dependency. Viscous strain rate hardening simply delays its appearance.

In the last years, interest has been shown in the combination of viscous material response and multiphysics couplings as means of regularization (see [5,8,38–48]). The target in these approaches is again to achieve a localization width, as a mean of regularization, that does not depend on the discretization, but rather on the material properties.

More specifically, in the case where the multiphysics processes are of diffusive nature (e.g. temperature diffusion, pore fluid diffusion, chemical decomposition and solute diffusion, dislocation diffusion), it is thought that the diffusion of the physical quantity linked to the mechanical softening via the multiphysics couplings, will also lead to a diffused strain localization profile in the material. In essence, the expectation is that the incorporation of more physics into the problem will remove the singularities in the strain localization profile predicted by the mechanics alone.

In the present work, the regularizing potential of an advection-diffusion multiphysics process is investigated. The aim is to determine under which conditions the regularization of the problem is achieved due to multiphysics and what the characteristics of the localization pattern are. In essence, we investigate whether the characteristic length introduced by the multiphysics couplings (see for example diffusion processes) can regularize the problem or not.

More precisely, we investigate the case of a single physics mechanism coupled to the mechanical response, in the simplified 1D setting, in order to draw meaningful conclusions considering strain localization and the regularization potential of the coupling terms. Even though there is one physics mechanism involved, the nature of the couplings in this mechanism is general. Thus, the physics process discussed in the paper may correspond to different multiphysics processes depending on the nature of the couplings involved. We mainly focus on the study of each coupling term separately in order to provide bounds in the coupling parameters, under which strain regularization is successful.

The paper is structured as follows: In Section 2, we present the material response under shear due to the combination of normal and shear stresses. The material elasto-viscoplastic mechanical response is coupled to the general form of an advection-diffusion equation (e.g. due to the percolation of pore fluid through the pores of a geomaterial). We take into account the effect of the different terms that are present in the final dispersion relation. In Section 3, we investigate the pure mechanical behavior (i.e. the mechanical strain softening and strain rate hardening). In Section 4, we investigate the influence of the multiphysics couplings on an elasto-plastic Cauchy continuum. In Section 5 we study the influence of multiphysical couplings in the presence of viscosity. We focus on the role of different terms in the advection-diffusion equation, namely: a) the diffusion coefficient, b) the total strain, c) the total strain rate, d) the effect of the amount of the diffusion physical quantity $\theta$ on the diffusion coefficient and on itself (coefficients A, B on the main text respectively). The effects of inertia are also taken into account. Next, in Section 7, we present a 1D numerical implementation of the above theoretical results by considering a fault gouge under shear taking the influence of the hydro-mechanical couplings into consideration. Finally, we discuss the theoretical and numerical results considering their implications when modeling failure under the effect of multiphysics couplings.

## 2. Governing equations

In this paper our goal is to understand under which conditions strain regularization and mesh independence can be achieved in the presence of a single multiphysics process. Thus, in order to keep the presentation short, we study here the simplest case considering the response of a 1D physics coupled mechanical system. Moreover, strain localization in two and three dimensions develops in the form of bands whose kinematics and equilibrium consideration do not differ significantly from the one dimensional case. The general equations

<figure><img src="image_4.png" /><figcaption>Fig. 1. Shearing of a 1D layer invariant in the $x_1$ direction. The layer is compressed by the application of a constant normal stress $\sigma_{22}$. Inside the layer a physical mechanism acts through advection and diffusion of a physical quantity affecting the stress distribution.</figcaption></figure>

of this system consist of the two way coupled stress balance and the advection-diffusion equations.

We use linear stability analysis to examine the stability of the uniform deformation equilibrium condition, according to the first method of Lyapunov (see [49]). In practice, starting from a homogeneous solution to the above system, we apply a small perturbation to it and we examine the evolution of its amplitude.

### 2.1. Stress balance

We describe the boundary value problem shown in Fig. 1, using the balance equations for the total stresses $\sigma_{ij}$

$$ \sigma_{ij,j} = \rho \ddot{u}_i, \quad i = 1, 2, 3 \text{ and } j = 1, 2, 3, \qquad (1) $$

where Einstein notation is used, $\rho$ is the density, $u_i$ is the displacement, $(\cdot)_j$ defines the spatial derivative with respect to $x_j$ and $(\cdot)$ denotes the second order time derivative.

The kinematic description of the problem is then completed by the following consideration for the total strains $\epsilon_{ij}$,

$$ \epsilon_{ij} = \frac{1}{2} (u_{i,j} + u_{j,i}). \qquad (2) $$

The total stress is assumed to consist of a purely mechanical part and a part resulting from the multiphysics coupling

$$ \sigma_{ij} = \sigma'_{ij} + \alpha_{ij} \theta, \qquad (3) $$

where $\theta$ is a quantity obeying an advection-diffusion equation, such as for example temperature variation or pore fluid pressure. The coupling between the total stresses and the physical quantity is achieved through the use of a tensor of coupling coefficients $\alpha_{ij}$. For the sake of simplicity, in the rest of our analysis, we will consider the coupling coefficients $\alpha_{ij}$ to be constant and independent of the physical quantity $\theta$ and the displacements $u_i$ to any order.

Taking into account that geometry and loading are uniform in two directions, the displacement fields $u_1, u_2$ are considered to vary only in the vertical direction $x_2$ (invariance in $x_1$). The system then reduces to

$$ \begin{aligned} \sigma_{12,2} &= \rho \ddot{u}_1, \\ \sigma_{22,2} &= \rho \ddot{u}_2, \end{aligned} \qquad (4) $$

where $u_1 = u_1(x_2), u_2 = u_2(x_2), \sigma_{12} = \sigma_{12}(x_2), \sigma_{22} = \sigma_{22}(x_2)$. In the next section, we complete the mathematical description of the boundary value problem, by introducing the constitutive law relating the rates of the effective stresses and strains.

### 2.2. Constitutive description

We consider a Perzyna elasto-viscoplastic constitutive model for the mechanical effective stress. This model describes the intrinsic rheological behavior of the material under creep and relaxation conditions
(see [34,36]). According to [50], the elasto-viscoplastic potential can
the be given by

$$
\Omega(F, \dot{\lambda}) = F - F^{-1}(\eta F_0 \dot{\lambda}), \qquad (5)
$$

where, η is the fluidity parameter [s⁻¹], and ˙λ = (⟨F⟩/ηF₀)ⁿ is the elasto-
viscoplastic multiplier. The quantity ⟨F⟩ in the Macaulay brackets is
equal to F when F ≥ 0 and 0 otherwise, where F is the yield surface.
Finally, F₀ denotes the initial value of the yield criterion at zero plastic
deformation. In the following example, we consider the yield surface
to be dependent on the effective values of the shear stress σ'₁₂,

$$
F = \sigma'_{12} - \tau_0 - H \epsilon_{12}^{vp}, \quad (6)
$$

where τ₀ is the initial yield strength of the material and H is the hard-
ening (softening) parameter of the material during the accumulation
of the visco-plastic strain ε¹². For a general viscoplasticity formulation
in the full stress space, we refer to [37]. Taking into account small
deflections, we consider an additive decomposition of the strain
increments into elastic and viscoplastic parts

$$
\dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^{(el)} + \dot{\epsilon}_{ij}^{(vp)}. \quad (7)
$$

We consider the following relationship between stresses and elastic
strains

$$
\sigma'_{ij} = M_{ijkl}^{(el)} \epsilon_{kl}^{(el)}, \quad (8)
$$

which, for $\sigma'_{12}$ becomes

$$
\sigma'_{12} = 2G\epsilon_{12}^{(el)}, \qquad (9)
$$

where, G is the shear modulus.

Given the form of the yield function, ε'_{12}^{vp} is the only viscoplastic strain, with the other strain components having only an elastic part. The viscoplastic strain rate is

$$
\epsilon_{12}^{vp} = \lambda \frac{\partial F}{\partial \sigma'_{12}}, \qquad (10)
$$

and in accordance to the Perzyna model for $n = 1$, we obtain

$$
\epsilon_{12}^{vp} = \frac{\langle F \rangle}{\eta F_0} \frac{\partial F}{\partial \sigma'_{12}} = \frac{\sigma'_{12} - \tau_0 - H \epsilon_{12}^{vp}}{\eta F_0}. \quad (11)
$$

In order to proceed with the stability analysis, we need to formulate
the equilibrium equation of the system in terms of the principal un-
knowns, i.e the displacements uᵢ, the physical quantity θ and their rates.
Here, we follow the process described in [37,51]. When the system is
under viscoplastic loading, F > 0, combining (11) with (7), rearranging
the terms and solving for the effective stress σ'₁₂, we obtain

$$
\sigma'_{12} = \frac{1}{1+h} \left( 2Gh\epsilon_{12} + 2\eta G f_0 \dot{\epsilon}_{12} - \eta f_0 \dot{\sigma}'_{12} + \tau_0 \right), \quad (12)
$$

where we set $f_0 = F_0/2G$ and $h = H/2G$. Next, we combine (12) with (3) and (4) to obtain the equilibrium equation of an elasto-viscoplastic layer including inertia forces,

$$
\frac{1}{1+h} \left( 2Gh\epsilon_{12,2} + 2\eta G f_0 \dot{\epsilon}_{12,2} - \eta f_0 \dot{\sigma}_{12,2} \right) + \alpha_{12} \theta_{,2} = \rho \ddot{u}_1. \quad (13)
$$

In addition, we make use of the first time derivative of (4) and
substitute displacements for the strains. Thus we obtain

$$
\begin{align}
& \frac{G}{1+h} (h u_{1,22} + \eta f_0 \dot{u}_{1,22}) + \left( \alpha_{12} + \frac{\eta f_0}{1+h} \dot{\alpha}_{12} \right) \dot{\theta}_{,2} + \frac{\eta f_0 \alpha_{12}}{1+h} \dot{\theta}_{,2} \nonumber \\
& \qquad = \rho \ddot{u}_1 + \frac{\eta f_0}{1+h} \rho \ddot{u}_1. \tag{14}
\end{align}
$$

We note that when η = 0, we obtain the wave equation for an
elastoplastic strain hardening medium. Moreover, we consider that
˙α₁₂ = 0, i.e. the coupling tensor is independent of the time. Linearizing
Eqs. (14) around the homogeneous solution (uᵗ, θᵗ) we obtain

$$
\bar{h}\bar{u}_{1,22} + \bar{\eta}\bar{u}_{1,22} + \bar{\alpha}_{12}\dot{\bar{\theta}}_{,2} + \bar{\eta}\bar{\alpha}_{12}\dot{\bar{\theta}}_{,2} = \bar{\rho}\bar{u}_1 + \bar{\eta}\bar{\rho}\bar{u}_1, \quad (15)
$$

where

$$
\begin{align}
\bar{\eta} &= \frac{\eta f_0}{1+h}, \\
\bar{h} &= \frac{h}{1+h}, \\
\bar{\alpha} &= \frac{\alpha_{12}}{G}, \\
\bar{\rho} &= \frac{\rho}{G} = \frac{1}{v_s^2}. \tag{16}
\end{align}
$$

and $\bar{u}_1$ is the perturbation of the displacement.

2.3. Advection-diffusion equation

The general form of the system of advection-diffusion equations
depending on the physical quantities θᵢ is given by

$$
(D^{(i)} (u_{m,n}, \theta_k) \theta_{i,j})_j - (v_j \theta_i), j + R^{(i)} (u_{m,n}, \dot{u}_{m,n}, \theta_k) = \dot{\theta}_i, \quad (17)
$$

where D⁽ᵏ⁾ is the diffusion coefficient of the quantity θᵢ, that can depend
on θᵢ and uᵢⱼ, R⁽ᵏ⁾ is the source (or sink) term of the same quantity and
vⱼ is the velocity field controlling the advection. Moreover, we consider
that the source (sink) term depends also on uᵢⱼ, ˙uᵢⱼ, θᵢ.

Linearizing the diffusion and source terms D⁽ᵏ⁾, R⁽ᵏ⁾ of the balance Eq. (17), using a Taylor series around the homogeneous solution vector (uᵗ, θᵗ), leads to the formulation

$$
\begin{align}
& D^{(i)} \bar{\theta}_{i,j,j} + \nonumber \\
& \frac{\partial D^{(i)}}{\partial u_{m,n}} \theta_{i,j}^* \bar{u}_{m,n,j} + \frac{\partial D^{(i)}}{\partial u_{m,n}} u_{m,n,j}^* \bar{\theta}_{i,j} + \nonumber \\
& \frac{\partial D^{(i)}}{\partial \theta_k} \theta_{i,j}^* \bar{\theta}_{k,j} + \frac{\partial D^{(i)}}{\partial \theta_k} \theta_{k,j}^* \bar{\theta}_{i,j} - \nonumber \\
& \qquad v_{j,j} \bar{\theta}_i - v_j \bar{\theta}_{i,j} + \nonumber \\
& \frac{\partial R^{(i)}}{\partial u_{m,n}} \bar{u}_{m,n} + \frac{\partial R^{(i)}}{\partial \dot{u}_{m,n}} \dot{u}_{m,n} + \frac{\partial R^{(i)}}{\partial \theta_k} \bar{\theta}_k = \dot{\theta}_i, \tag{18}
\end{align}
$$

where the dependencies of the functions D⁽ᵏ⁾ (uₘ,ₙ, θₖ) and R⁽ᵏ⁾ (uₘ,ₙ, ˙uₘ,ₙ, θₖ) have been omitted for the sake of brevity. Collecting the terms with the same perturbation components, limiting our analysis to one physical process and to one dimension simplifies (18) to

$$
D \bar{\theta}_{22} + A \bar{\theta}_{2} + B \bar{\theta} + U \bar{u}_{1,22} + V \bar{u}_{1,2} + W \bar{u}_{1,2} = \dot{\bar{\theta}}, \quad (19)
$$

where

$$
\begin{align}
A &= \left( \frac{\partial D}{\partial u_{1,2}} u_{1,22}^* + 2 \frac{\partial D}{\partial \theta} \theta_{,2}^* - v_2 \right), \nonumber \\
B &= \left( \frac{\partial R}{\partial \theta} - v_{2,2} \right), \nonumber \\
U &= \frac{\partial D}{\partial u_{1,2}} \theta_{,2}^{*}, \quad V = \frac{\partial R}{\partial u_{1,2}}, \quad W = \frac{\partial R}{\partial \dot{u}_{1,2}}. \tag{20}
\end{align}
$$

Furthermore, the term u<sub>1,22</sub><sup>*</sup> is zero because we wish to characterize the
stability of uniform deformations.

The existence of the unique solution of the system of linearized
equations comprising of Eqs. (15) and (19) that satisfy the boundary
conditions, is given by (¯u₁, ¯θ),

$$
\begin{bmatrix} \bar{u}_1 \\ \bar{\theta} \end{bmatrix} = \begin{bmatrix} g \\ y \end{bmatrix} e^{st+ikx_2}, \qquad (21)
$$

where s = s<sub>r</sub> + s<sub>i</sub> is the Lyapunov coefficient, i = √-1. Its real
part, s<sub>r</sub> controls the growth of the perturbation w.r.t. time. When an
imaginary part s<sub>i</sub> is also present, then a traveling pulse will occur
(see [37] and Appendix B.1). Considering that k is the wavenumber
corresponding to the wavelength of the solution λ = 2π/k, then we obtain
its traveling velocity c = s<sub>i</sub>/k. Introducing the solution to the linearized
balance equations Eqs. (15), results in the following system involving
the solution constants,

$$
\begin{bmatrix} \bar{\eta} \bar{\rho} s^3 + \bar{\rho} s^2 + \bar{\eta} k^2 s + \bar{h} k^2 & -ik\bar{\alpha} (\bar{\eta} s + 1) \\ U k^2 - \text{Vik} - \text{Wiks} & s + Dk^2 - \text{Aik} - B \end{bmatrix} \begin{bmatrix} g \\ y \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}. \quad (22)
$$
<figure><img src="image_5.png" /><figcaption>Fig. 2. The roots of the characteristic polynomial (Lyapunov coefficients, s) in a Cauchy medium: (a) The real positive roots for: (i) plain Cauchy (s₁), (ii) rate dependence no inertia (s⁺), (iii) rate dependence and inertia (sᵢ). Cases (i,iii) are mesh dependent. (b) The remaining roots of case (iii), (sᵢᵢ, sᵢᵢᵢ) become complex conjugate when λ₀ is reached. (c) The critical wavenumber k₀ = 1/λ₀ depends on the softening parameter ĥ. (d) When multiphysics couplings are considered with (ii), the conjugate roots, sᵢᵢ, sᵢᵢᵢ (black line), become positive near λ = 0 and present a global maximum at λ = 0, over the real root (sᵢᵢᵢ, red curve). Mesh dependence occurs.</figcaption></figure>

Apart from the trivial solution [g, y]ᵀ = [0, 0], the system admits an infinity of solutions when its determinant becomes zero,

$$
\begin{align}
\Delta &= \bar{\eta}\bar{\rho}s^4 + (D\bar{\eta}\bar{\rho}k^2 - iA\bar{\eta}\bar{\rho}k - B\bar{\eta}\bar{\rho} + \bar{\rho})s^3 \nonumber \\
&\quad + [(D\bar{\rho} + \bar{\eta} + \bar{\alpha}\bar{\eta}W)k^2 - iA\bar{\rho}k - B\bar{\rho}]s^2 \nonumber \\
&\quad + [D\bar{\eta}k^4 + i\bar{\eta}(\bar{\alpha}U - A)k^3 + (\bar{\alpha}\eta V + \bar{\alpha}W - B\bar{\eta} + \bar{h})k^2]s \nonumber \\
&\quad + D\bar{h}k^4 + i(\bar{\alpha}U - A\bar{h})k^3 + (\bar{\alpha}V - B\bar{h})k^2 \nonumber \\
&= 0. \tag{23}
\end{align}
$$

The stability of the system of equations (i.e. whether the solution converges, diverges or is indeterminate w.r.t. the homogeneous solution) is controlled by the sign of the real part of s. Due to the very large number of parameters, it is difficult to analyze the stability of the system as a whole. Therefore, in our analysis, we will discuss a combination of special cases taking into account the effects of individual terms together with appropriate examples.

### 3. Uncoupled stress balance

We will start our analysis by first considering the case of the uncoupled mechanical system (ā₁₂ = 0). In this case the determinant reduces to a third order polynomial

$$ \Delta = \bar{\eta}\bar{\rho}s^3 + \bar{\rho}s^2 + \bar{\eta}k^2s + \bar{h}k^2, \quad (24) $$

whose coefficients are all positive with the exception of the constant term (ħk²). The discriminant of the cubic polynomial of (24) is equal to

$$ D_v = k^2 \bar{\rho} (\bar{\eta}k^2(\bar{\eta}\bar{\rho} + 18\bar{h}\bar{\rho}\bar{\eta}) - 4\bar{h}\bar{\rho}^2 - 4\bar{\eta}^4 k^4). \quad (25) $$

We note that the sign of (25) depends on the second term of the product. The second term is a quartic polynomial in k whose leading coefficient is negative. Therefore, for large values of the wavenumber k the value of the discriminant is negative. This indicates that the polynomial has one real root and two conjugate complex roots.

#### 3.1. Positive root

For large wavenumbers (k → ∞), based on the rule of signs, we establish the existence of one positive real root, which is equal to

$$ s_0 = -\frac{h}{f_0 \eta}. \quad (26) $$

We present this result in Fig. 2(a) and compare with the simple cases of a rate dependent medium in the absence of inertia (see [32]) and the classical rate independent medium in the presence of inertia. We note that in the absence of inertia, the value of the positive root is independent of k. Moreover, we establish that in the classical rate independent case (presence of inertia and absence of viscosity), s goes to infinity (s → ∞) when the wavelength tends to zero (λ → 0).

This indicates that without the consideration of inertia any initial defects will be propagated without a change in their width, while with inertia strain localization on a mathematical plane and mesh dependence of the numerical results will take place.

We note that even though the introduction of viscosity in the equilibrium equation considering inertia leads to a characteristic length (see [36,52]), this is not enough for the regularization of strain and mesh independence. This was also shown to be the case in [35,37,42].

#### 3.2. Remaining roots

For high values of the wavenumber k, we turn our attention to the sign of the real part of s in the complex remaining roots of (24). We decompose s into its real and imaginary parts and substitute into the third root polynomial of (24). The solution for the real part of the complex roots is a cubic polynomial

$$ -8\bar{\eta}\bar{\rho}s_r^3 - 8\bar{\rho}s_r^2 - \left(\frac{2\bar{\rho}}{\bar{\eta}} + 2\bar{\eta}k^2\right)s_r - k^2(1 - \bar{h}) = 0. \quad (27) $$

Making use of the Routh-Hurwitz criterion (see [53,54]), it results that the above can only have negative solutions, so that any complex roots must have negative real parts. Dividing all terms in (27) by k² and taking the limit k → ∞ yields

$$ s_r = -\frac{1}{2\eta f_0}. \quad (28) $$
Thus, at large wavenumbers, we note that the complex conjugate solutions have negative real parts and the corresponding perturbations decay with time. When we decrease the value of the wavenumber, we reach a value, $k_0$, for which the discriminant changes sign and becomes positive. In this case, according to the rule of signs, the roots of the polynomial become real and negative, leading the corresponding perturbations to decay. The value for which the discriminant changes sign is given by

$$k_0 = \sqrt{\frac{18\bar{h}\bar{\rho} + \bar{\rho} + \sqrt{(18\bar{h} + 1)^2 - 64\bar{h}}}{8\bar{h}^2}} \quad (29)$$

In the left part of Fig. 2(c), we present the critical wavelength for which the discriminant changes sign from negative to positive w.r.t the softening variable $\bar{h}$ for the specific choice of mechanical parameters $\bar{\rho} = 0.01$, $\bar{\eta} = 1.0$. In the right part of Fig. 2(c), we present the remaining roots of the system, which transition from complex conjugate to real values. In accordance to the main theoretical result, their real part is always negative.

From the above analysis, we show that under the influence of inertia and viscosity the solution to the system of equations diverges from the homogeneous solution (fixed point). Moreover, the part of the solution that diverges the fastest is that of the zero wavelength. Therefore, the solution localizes on a mathematical plane (i.e. a zero thickness shear band is formed) and mesh dependence will occur. We conclude that inertia and viscosity by themselves are not capable of providing an internal length that regularizes the problem of strain localization. For this reason we will consider in the next sections the coupling of the mechanical equation with the advection-diffusion process resulting in the general system of (22).

## 4. Multiphysics coupling with rate-independent Cauchy continuum

In this section we will take into account the two-way multiphysics coupling imposed in a Cauchy continuum with a rate-independent elastoplastic mechanical behavior and the corresponding advection-diffusion equation. In Section 4.1 we will first investigate the case of multiphysics couplings in the absence of inertia. Next, in Section 4.2 we will investigate the result in the presence of inertia.

### 4.1. Cauchy rate independent medium in the absence of inertia

In the absence of inertia, the system of (22) simplifies to

$$\begin{bmatrix} \bar{h}k^2 & -ik\bar{\alpha} \\ Uk^2 - \text{Vik} - \text{Wiks} & s + Dk^2 - \text{Aik} - B \end{bmatrix} \begin{bmatrix} g \\ y \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}. \quad (30)$$

We note that in the system of (30), the Lyapunov coefficient is found only in the advection-diffusion equation. The mechanical part of the system is time insensitive. Because of this, we cannot assert that the system's homogeneous solution is stable. However, we can show under which conditions the homogeneous solution is unstable (i.e. the Lyapunov coefficient has positive real part $s_r$). The system determinant is equal to

$$\Delta = [(\bar{\alpha}W + \bar{h})k^2] s + D\bar{h}k^4 + i(\bar{\alpha}U - A\bar{h})k^3 + (\bar{\alpha}V - B\bar{h})k^2, \quad (31)$$

meaning that the determinant of the system becomes zero when

$$s = - \frac{D\bar{h}k^2 + i(\bar{\alpha}U - A\bar{h})k + (\bar{\alpha}V - B\bar{h})}{(\bar{\alpha}W + \bar{h})}. \quad (32)$$

In order to decide whether the solution of zero wavelength will dominate (strain localization on a mathematical plane) the rest of the solutions of system (30), we need to find the value of the real part of $s$ as the wavelength, $\lambda$, tends to zero (respectively the wavenumber, $k$, tends to infinity). In the case of small wavelengths (large wavenumbers $k \to \infty$) we obtain

$$s = - \frac{D\bar{h}}{\bar{\alpha}W + \bar{h}} \lim_{\lambda \to 0} \frac{1}{\lambda^2}. \quad (33)$$

**Table 1**

Synopsis of the theoretical results for the stability of the homogeneous solution (h.s) in a multiphysics coupled rate independent Cauchy continuum without inertia. We present the limit of the real part of the Lyapunov coefficient, $s$, for small wavelengths ($\lambda \to 0$). We focus on the cases of the mechanical hardening (respectively softening) $\bar{h} > 0$, $(\bar{h} < 0)$ and their combination with the sign and magnitude of the deformation rate sensitivity term $\bar{\alpha}W$. $\bar{\alpha}W > 0$, (respectively $\bar{\alpha}W < 0$) leads to increase (decrease) of the physical quantity $\theta$ for different strain rates). Based on the first method of Lyapunov, we can only assert strain localization on a mathematical plane when $s \to +\infty$.

<table><thead><tr><th>h&#773;</th><th>ᾱW</th><th>s</th><th>Comment</th></tr></thead><tbody><tr><td>&gt;0</td><td>&gt; -h&#773;</td><td>-∞</td><td>h.s of advection-diffusion stable, no decision on strain localization</td></tr><tr><td>&gt;0</td><td>&lt; -h&#773;</td><td>+∞</td><td>h.s of advection-diffusion unstable, strain localization on a mathematical plane</td></tr><tr><td>&lt;0</td><td>&gt; -h&#773;</td><td>+∞</td><td>h.s of advection-diffusion unstable, strain localization on a mathematical plane</td></tr><tr><td>&lt;0</td><td>&lt; -h&#773;</td><td>-∞</td><td>h.s of advection-diffusion stable, no decision on strain localization</td></tr></tbody></table>

Considering that $D > 0$, we notice that strain localization on a mathematical plane will occur for specific values of the normalized hardening (respectively softening) modulus, $\bar{h} > 0$, $(\bar{h} < 0)$, and the sensitivity of the source term $R$ to the deformation rate $\bar{\alpha}W$.

More specifically, considering the quantity $\bar{\alpha}W$, we note that when the term $\bar{\alpha}W$ is positive the rate of the physical quantity $\theta$ in the diffusion equation increases for faster shearing velocities, else it decreases. Thus, we define that when $\bar{\alpha}W > 0$ (respectively $\bar{\alpha}W < 0$), there is a positive (negative) feedback from the mechanical system to the advection-diffusion equation.

Next, we summarize in Table 1, the sign of $s$, w.r.t. the relative difference, $\bar{h} - \bar{\alpha}W$. We note that for small wavelengths ($\lambda \to 0$) and $s \to +\infty$, the solution of the system (30) localizes on a mathematical plane. In the case where $s \to -\infty$, the fastest growing solution component is given by the limit of $s$ for large wavelengths $\lambda \to \infty$, and no strain localization on a mathematical plane will take place.

We conclude from (33) and Table 1 that strain localization on a mathematical plane and mesh dependence of the numerical results occur when:

* Mechanical hardening, $\bar{h} > 0$, is combined with strong negative feedback $\bar{\alpha}W < -\bar{h}$.

* Mechanical softening $\bar{h} < 0$, is combined with weak negative or positive feedback $\bar{\alpha}W > -\bar{h}$.

Here, we have assumed that the diffusion coefficient is positive. For negative diffusion the signs in the table are inverted. In the absence of diffusion the real part of the Lyapunov coefficient, $s$, is constant and the imaginary part of $s$ grows with increasing wavenumber $k$.

### 4.2. Cauchy rate independent medium in the presence of inertia

When inertia is taken into account, we obtain the full form of the determinant

$$\Delta = \bar{\rho}s^3 + \bar{\rho}[Dk^2 - iAk - B]s^2 + [(\bar{\alpha}W + \bar{h})k^2]s + D\bar{h}k^4 + i(\bar{\alpha}U - A\bar{h})k^3 + (\bar{\alpha}V - B\bar{h})k^2. \quad (34)$$

For large values of the wavenumber, $k \to \infty$, ($\lambda \to 0$), the determinant can be approximated as

$$\Delta = \bar{\rho}s^3 + \bar{\rho}Dk^2s^2 + (\bar{\alpha}W + \bar{h})k^2s + D\bar{h}k^4. \quad (35)$$

In this case the imaginary part in the second order term, $s^2$, and the constant term do not contribute in the polynomial. Moreover, we assume that mechanical softening takes place, $\bar{h} < 0$. Next, we use the rule of signs, to obtain the presence of exactly one root with a positive real part for large values of the wavenumber ($k \to \infty$). Isolating the highest power terms of $s, k$, we note that when $k \to \infty$, the equality
can only hold when $s \to \infty$. The values of the three roots at small wavelength $\lambda \to 0$, (large wavenumbers, $k \to \infty$) tend to $\lim_{k \to 0} s_i = \lim_{k \to \infty} s_i = (-1)^i \infty$, $i = 1, ..., 3$.

When small values of the wavenumber are taken into account, $k \to 0$, the determinant becomes

$$ \Delta = \bar{\rho}s^2(s - B), \qquad (36) $$

leading to

$$ \lim_{k \to 0} s_{1,2} = \lim_{\lambda \to \infty} s_{1,2} = 0 \quad \text{and} \quad \lim_{k \to 0} s_3 = \lim_{\lambda \to \infty} s_3 = B, \qquad (37) $$

where $B = \frac{\partial R}{\partial \theta} - v_{2,2}$. This indicates that the third root, $s_3$, is negative at $k \to \infty$. Moreover, it does not tend to zero at $k \to 0$, but it is either a positive or negative finite number.

In Fig. 2(d), we present an application of the above theoretical results in the special case where $\bar{h} = -0.25$, $\rho = 0.001$, $D = \bar{\alpha} = B = W = 1$, $A = U = V = 0$. The non zero terms correspond to the coupling term in the stress balance equation ($\bar{\alpha}$), and the dependence of the source term in the advection-diffusion equation on $\theta$ ($B$) and on the strain rate ($W$). The real root is shown (red curve) along with the real part of the complex conjugate roots (black curve).

We notice that a global maximum, $s \to +\infty$ is present as $\lambda$ goes to zero leading to strain localization on a mathematical plane and mesh dependence of the numerical results. By superposition of the 3 roots, we note here the existence of a local maximum for an intermediate value of $\lambda$. In such a case, when a numerical analysis with finite elements takes place, the strain profile will show signs of localizing on a width equal to the wavelength of the local maximum when the size of the elements used is larger than that. For finer discretizations, mesh dependence of the numerical results will occur (localization on a mathematical plane).

## 5. Multiphysics coupling with viscosity

In this section we will discuss the effect of multiphysics couplings and viscosity on the stability of the homogeneous solution during shearing of the material. We will show that, depending on the values of the different coupling parameters in the diffusion equation, the consideration of viscosity and multiphysics can lead to a) a stable homogeneous solution (strain regularization), b) strain localization on a finite wavelength (unstable homogeneous solution, strain regularization) and c) strain localization on a mathematical plane. Moreover, we explain the implications of our mathematical analysis for the expected numerical results.

Due to the separate influence of the different terms present in the coupled diffusion equation describing the evolution of the multiphysics coupling (see (18)), our analysis has multiple practical applications covering a broad spectrum of interactions between the response of the material and multiphysics coupling. In each subsection below, we provide an application, that is an example of a material whose behavior falls under the examined coupling description. Finally, we consider in Section 5.1.3 a particular numerical application to one of the cases discussed here.

### 5.1. In the absence of inertia

In the absence of inertia, the determinant simplifies to

$$ \Delta = k^2(Rs^2 + P(k)s + Q(k)), \qquad (38) $$

where

$$ \begin{aligned} R &= \bar{\eta}(1 + \bar{\alpha}W), \\ P(k) &= [D\bar{\eta}k^2 + i\bar{\eta}(\bar{\alpha}U - A)k + (\bar{\alpha}\eta V + \bar{\alpha}W - B\bar{\eta} + \bar{h})], \\ Q(k) &= (D\bar{h}k^2 + i(\bar{\alpha}U - A\bar{h})k + (\bar{\alpha}V - B\bar{h})). \end{aligned} \qquad (39) $$

We note that $k^2$ is a common factor in the above dispersion relation. The remaining factor is then a quadratic polynomial in $s$, with roots

$$ s(k) = \frac{-P(k) \pm \sqrt{P(k)^2 - 4RQ(k)}}{2R}. \qquad (40) $$

These roots control the stability of the system.

In order to study the localization width in this case, we focus on the behavior of the solution $s$ for large and small values of the wavenumber $k$. Taking the limit of the solution for large values of the wavenumber, $k \to \infty$, results in an indeterminate form, $(\infty \pm \infty)$, in the numerator. Rearranging the terms in the fraction to resolve the indeterminacy, and further noticing that $\lim_{k \to \infty} P(k) > 0$, considering $D > 0$, $\bar{\eta} > 0$, we obtain

$$ \lim_{\lambda \to 0} s_{1,2} = \lim_{k \to \infty} \frac{P(k)}{2R} \left( -1 \pm \sqrt{1 - \frac{4RQ(k)}{P^2(k)}} \right). \qquad (41) $$

Finally, for the different signs this limit results in

$$ \lim_{\lambda \to 0} s_1 = \lim_{k \to \infty} s_1 = -\frac{\bar{h}}{\bar{\eta}}, \qquad (42) $$

$$ \lim_{\lambda \to 0} s_2 = \lim_{k \to \infty} s_2 = -\lim_{k \to \infty} \frac{P(k)}{R}, \qquad (43) $$

with the first being positive for softening material and the sign of the second depending on the sign of the term $R$.

In the cases where $\lim_{\lambda \to 0} s_1 = -\frac{\bar{h}}{\bar{\eta}}$ the Lyapunov coefficient obtains a finite positive value at zero wavelength. Moreover, this is a local maximum of the dispersion relation and does not imply strain localization on a mathematical plane. In fact, strain localization on a mathematical plane can be avoided if we can find larger values of the Lyapunov coefficient, $s$, for different values of the wavelength $\lambda \neq 0$.

For this reason, we compare with the values of the Lyapunov coefficient for large wavelengths, $\lim_{\lambda \to \infty} s_{1,2}$, and check if it is greater than $\lim_{\lambda \to 0} s_1$. In the case of large wavelengths $\lambda \to \infty$ the behavior of the solution depends on the quadratic polynomial term of Eq. (38),

$$ \lim_{\lambda \to \infty} s = \lim_{k \to 0} s_{1,2} = \frac{-P(0) \pm \sqrt{P(0)^2 - 4RQ(0)}}{2R}. \qquad (44) $$

If $\lim_{\lambda \to \infty} s_{1,2} > \lim_{\lambda \to 0} s_1$ then no strain localization on a mathematical plane will take place.

When we take into account the second root, $s_2$, of the dispersion relation in (43), we establish that the value of the Lyapunov coefficient for small wavenumbers, depends on the relationships between $R$, $P(0)$, $Q(0)$, which are dependent on the mechanical and physical equations of the system and their couplings i.e the system parameters $\bar{\eta}$, $D$, $\bar{\alpha}W$, $\bar{\alpha}U$, $A$, $B$, $U$, $V$.

In the rest of the section, in order to decide when strain localization on a mathematical plane is possible, we examine the separate influence of the different terms in Eqs. (43), (44) and compare between them. Throughout our analysis, we will consider that the softening modulus lies in the range $0 > \bar{h} > -1$, and that the diffusion and viscosity coefficients are always positive, $D > 0$, $\bar{\eta} > 0$. In Figs. 3(a) to 3(d), we present the Lyapunov coefficients for different values of the system parameters. Finally, we summarize our findings in Table 2.

### 5.1.1. Dependence of the diffusion coefficient on the strain ($\bar{\alpha}U, \bar{\alpha}V = 0, \bar{\alpha}W = 0, A = 0, B = 0$)

The term $U$ stands for the dependence of the diffusion coefficient on the strain and is coupled to $\theta$. In practice we consider a porous material whose diffusion coefficient increases as the volume of the voids, which are indicative of the permeability, expands. As an example, this phenomenon is present in materials exhibiting hydraulic and chemical couplings (see also [26,55,56]).

If the term $U$ is present, the term $A$, which is indicative of the advection, is likely to be as well, see (18) and the definitions in (20), but since more parameters are included in the latter, we consider them separately.
<figure><img src="image_6.png" /><figcaption>Fig. 3. First row: Real part of the Lyapunov coefficient $Re(s)$ w.r.t. the wavelength $\lambda$. Second row: Imaginary part of the Lyapunov coefficient $Im(s)$ w.r.t. the wavelength $\lambda$. In (a,b,c,d,e) there exist values for each parameter such that the maximum $Re(s)$ lies at $\lambda \neq 0$, strain regularization and mesh independence develops. In (a,d), $Im(s) \neq 0$ when $maxRe(s)$ lies at intermediate $\lambda$ and a traveling strain localization of finite thickness occurs. In (b,c) for specific selection of the parameters $\bar{a}V$, $\bar{a}W$, $maxRe(s)$ lies at $\lambda = 0$ where $Im(s) = 0$, stationary strain localization on a mathematical plane takes place and mesh dependence.</figcaption></figure>

When $R > 0$, $\bar{\eta}D > 0$, $\bar{a}W = 0$, the real part of the positive root shows a local maximum $\lim_{\lambda \to 0} s = -\frac{\bar{h}}{\bar{\eta}}$. We note that for large wavelengths, $\lambda \to \infty$, a larger value of the Lyapunov coefficient is possible for all values of $\bar{a}U$, see Fig. 3(a). Moreover in Fig. 3(a), we establish that a global maximum exists for a value of the wavelength $\lambda$ different than zero. This means that strain localization on a mathematical plane is no longer possible and strain will localize on a plane of finite thickness that is independent of the discretization we impose in the numerical solution.

In Fig. 3(a), we show that an imaginary part is also present and exhibits a local maximum. The sign of the term $\bar{a}U$ does not affect the magnitude of the real or imaginary part of the solution, but changes the sign of its imaginary part. The presence of the imaginary part will lead to the appearance of a traveling strain localization profile inside the domain of the problem in the numerical analyses with finite elements (see also [37,42]).

### 5.1.2. Dependence of the source term in the advection-diffusion equation on the strain ($\bar{a}V$, $\bar{a}U = 0$, $\bar{a}W = 0$, $A = 0$, $B = 0$)

The term $V$ stands for the dependence of the source term in the advection-diffusion equation on the strain. One practical example of a multi-physics coupling that has been proposed to be affected by the level of strain inside the material, is the mechanism of dynamic strain aging (see [57]). In this mechanism in particular, the diffusion of the dislocation density inside the material is affected by the strain level. This leads to a particular characterization of the failure and strain localization patterns inside the material.

Considering such a material, the dispersion relation presents two roots are present, whose limits for small wavelengths are given in Table 2. For large wavelengths the limits are

$$ \lim_{\lambda \to \infty} s_{1,2} = - \frac{\bar{a}\bar{\eta}V + \bar{h} \pm \sqrt{(\bar{a}\bar{\eta}V + \bar{h})^2 - 4\bar{a}\bar{\eta}V}}{2\bar{\eta}}. \quad (45) $$

We present the values of the Lyapunov coefficient in Fig. 3(b). More specifically, in Fig. 3(b), we observe that for negative values of $\bar{a}V$, the large wavelength limit is a global maximum for the real part of the positive root. Thus strain localization on a mathematical plane and mesh dependence of the numerical results is not possible in this case. When $\bar{a}V > 0$, we find that the global maximum of the Lyapunov coefficient is found at $\lambda = 0$. Thus, strain localization on a mathematical plane and mesh dependence of the numerical results will take place.

In Fig. 3(b), we note that the Lyapunov coefficient has a non zero imaginary part only for positive values of $\bar{a}V$ when its real part is negative. In this case, the solutions of the Lyapunov coefficient are conjugates with the same real parts and opposite imaginary parts (see also Fig. 3(b)).

### 5.1.3. Dependence of the source term in the advection-diffusion equation on the strain rate ($\bar{a}W$, $\bar{a}U = 0$, $\bar{a}V = 0$, $A = 0$, $B = 0$)

The term $W$ stands for the dependence of the source term in the advection-diffusion equation on the strain rate. This is a case with particular importance for geomechanics applications such as the shearing of fluid saturated fault gouges (see also [19,58]), hydraulic fracturing (see [59]), fluid flushing (see [60,61]), earthquake control (see [62,63]) and it will be further highlighted with the help of a numerical example in Section 7.

Two roots are present, whose limits for small wavelengths are given in Table 2. In the case of small wavelengths $\lambda \to 0$, we note that the value of the Lyapunov coefficient for the second root depends on the sign of the source term $R = 1 + \bar{a}W$. In particular we get

$$ \lim_{\lambda \to 0} s_1 = -\frac{\bar{h}}{\bar{\eta}}, \quad (46) $$

$$ \lim_{\lambda \to 0} s_2 = -\text{sign}(1 + \bar{a}W) \infty. \quad (47) $$

Furthermore, for large wavelengths the solutions of the Lyapunov coefficient are

$$ \lim_{\lambda \to \infty} s_{1,2} = - \frac{\bar{a}W + \bar{h} \pm |\bar{a}W + \bar{h}|}{2\bar{\eta}(\bar{a}W + 1)}. \quad (48) $$

Finally, we note that the solutions are anti-symmetric on the real line yielding

$$ \lim_{\lambda \to \infty} s_1 = 0, \\ \lim_{\lambda \to \infty} s_2 = - \frac{\bar{a}W + \bar{h}}{\bar{\eta}(\bar{a}W + 1)}. \quad (49) $$

In this case we note that strain localization on a mathematical plane and mesh dependence of the numerical results depend on the values of $\bar{a}W$. In particular, we describe three localization regimes,

* $\bar{a}W < -1$: Global maximum value of the Lyapunov coefficient, $\lim_{\lambda \to 0} s_1 \to +\infty$ at $\lambda \to 0$. The strain will localize on a mathematical plane and mesh dependence of the numerical results will take place.

* $-1 < \bar{a}W < 0$: We note that $\lim_{\lambda \to 0} s_1 \le \lim_{\lambda \to \infty} s_2$, which results in a global maximum at $\lambda \to \infty$. Strain regularization will take place and no mesh dependence will be observed in the numerical results.

* $\bar{a}W > 0$: Global maximum value of the real part $\lim_{\lambda \to 0} s_1 = -\frac{\bar{h}}{\bar{\eta}}$ at $\lambda \to 0$. The strain will localize on a mathematical plane and mesh dependence of the numerical results will take place.
We present the above conclusions in Fig. 3(c), which shows the real and imaginary parts of the Lyapunov coefficient for different wavelength values $\lambda$.

### 5.1.4. Effect of advective flow and dependence of the diffusivity on the strain and $\theta$ ($A$, $\bar{a}U = 1$, $\bar{a}V = 0$, $\bar{a}W = 0$, $B = 0$)

The practical importance of this term has been discussed already in Section 5.1.1. Nevertheless, this term depends also on the physical quantity $\theta$ itself. This is a common occurrence in different scenarios in geophysics, such as thermal diffusion in the earth's crust and faults (see [64–66]) and diffusion of chemicals under chemical decomposition during an earthquake (see [26,67,68]).

The term $\bar{a}$ is necessary to coupling and its absence means that a constant root with a value of $-\frac{\hbar}{\bar{\eta}}$ is present. It will, therefore, always be assumed present. In addition, both $U$ and $A$ express the variation of the diffusion coefficient with the strain, so $\bar{a}U$ will also be assumed to be present and equal to 1. The limits for small and large values of the wavelength are not affected by the presence of $A$.

The presence of $\bar{a}U$, ($\bar{a}U \neq 0$) leads to the presence of critical points for intermediate values of the wavelength $\lambda > 0$ in the real part of the Lyapunov coefficient $s$. The number of critical points depends on the values of $\bar{a}U$ and $A$. This is also shown in Fig. 3(d). More specifically, when $\lambda > 0$, for large values of $\bar{a}U$ (e.g. $\bar{a}U = 1.0$) and low values of $A$ (e.g. $A = -1$) two critical points are present. The first one is an inflection point and constitutes a global minimum and the second one is a global maximum. For larger values of the parameter $A$ the inflection point disappears, while the global maximum increases in magnitude.

Next, we consider the imaginary part in Fig. 3(d). We note that the imaginary part of $s$ is affected more than the real parts during the parametric investigation of $A$. Furthermore, we note that the existence of an imaginary part for all values of the wavelength depends on the presence of the $\bar{a}U$ term. For large values of the parameter $\bar{a}U$, (e.g. $\bar{a}U = 1$) the response of the imaginary part depends on the values of the parameter $A$.

In particular low values of the parameter $A$ lead to one intermediate critical wavelength, where a discontinuity of the imaginary part of the Lyapunov coefficient takes place (see Fig. 3(d)). This critical wavelength in the imaginary part of the Lyapunov coefficient coincides with the position of the inflection point in the diagram of the real part of the Lyapunov coefficient (see Fig. 3(d)). For large values of the parameter $A$ the discontinuity of the imaginary part ceases and a global minimum takes its place. The positions of the global maxima and minima in the real and imaginary parts of the Lyapunov coefficient do not coincide (see Fig. 3(d)).

From the above analysis, we conclude that the global maximum of the real part of the Lyapunov coefficient lies at an intermediate value of the wavelength $\lambda$ and coincides with a non zero imaginary part. A traveling band of finite width will develop leading to strain regularization and no mesh dependence.

### 5.1.5. Effect of divergence of the advective flow and dependence of the source term for $\theta$ on itself ($B$, $\bar{a}U = 0$, $\bar{a}V = 0$, $\bar{a}W = -0.5$, $A = 0$)

Finally, we study the effect of the variation of the source term with the physical quantity $\theta$ and the compressible part of the advective flow by varying the parameter $B$. We consider the presence of the coupling term $\bar{a}$ and a source term $W$, (dependence on the strain rate, $\bar{a}W = -0.5$). This is an important application for problems involving a source and advective terms (e.g. hydraulic fracturing and hydraulic flushing in geomechanics [59–61], laser and frictional welding in metals [69,70] among others).

We find that $B$ has no effect on the limits for small wavelengths, but affects those for large wavelengths. More precisely we get

$$ \lim_{\lambda \to \infty} s_{1,2} = - \frac{\bar{a}W - B\bar{\eta} + \hbar \pm \sqrt{(\bar{a}W - B\bar{\eta} + \hbar)^2 + 4B\hbar(\bar{a}\bar{\eta}W + \bar{\eta})}}{2(\bar{a}\bar{\eta}W + \bar{\eta})} \quad (50) $$

From (50), we conclude that $B$ has an effect on the magnitude of the real part of the Lyapunov coefficient and can lead to the presence of an imaginary part. Moreover, it can change the sign of the real part of the root and, therefore, stability.

Contrary to the limit of large wavelengths, $\lambda \to \infty$, $k \to 0$, the limit for small values of the wavelength, $\lim_{\lambda \to 0} s_{1,2}$, remains positive for all values of $B$. Thus, regularization is only possible when the limit at $\lim_{\lambda \to \infty} s_{1,2} > \lim_{\lambda \to 0} s_{1,2}$. An example is shown in Fig. 3(e) for $\bar{a}W = -0.5$. In that case an imaginary part is not present.

### 5.1.6. Summary of the results

We summarize the main theoretical results of Section 5.1 in Table 2. More specifically, we present in each row of the Table the term we studied in each case. In the second and third columns of the Table, we present the results of the analysis in the cases of small $\lambda \to 0$ and large wavelengths $\lambda \to \infty$ respectively.

In each case we show the values for both roots of the characteristic polynomial. We note that in the cases of $\bar{a}V$ strain localization on a mathematical plane is always possible for negative values of $\bar{a}V < 0$. Moreover, in the cases of $\bar{a}W$, $B$ strain localization is always possible for strongly negative values of $\bar{a}W < -1$. For intermediate values of $\bar{a}W$, $-1 < \bar{a}W < 0$, strain regularization takes place because the Lyapunov coefficient for large wavelengths, $\lambda \to \infty$, is greater that that of small wavelengths $\lambda \to 0$. Finally, when $\bar{a}W > 0$ strain localization on a mathematical plane takes place.

## 5.2. In the presence of inertia

Up until now we have considered in our mathematical analysis cases of multiphysics couplings where the influence of the inertia is negligible. This is true in cases where the dimensions of the part of the material that localize are small, or when the loading velocity is small. This is generally the case for the material that lies away from the advancing front of the localization (see also [65]). In the case where we are interested in the near region of the advancing localization front, e.g. the tip of a crack or fault, the influence of inertia on the results of the mathematical analysis might become important for small times (see also [65,71]).

In what follows we will take into account the presence of inertia in our mathematical analysis. We note that in the presence of inertia the degree of the characteristic polynomial becomes large and a mathematical analysis of the determinant, its roots and their signs becomes complicated.

To begin with, we examine the limit of the real part of the Lyapunov coefficient, $s$, as the wavelength, $\lambda$, tends to infinity. In that limit two of the roots tend to zero, one to $B$ and one to $-1/\bar{\eta}$. Next we examine the limit of the real part of $s$ as $\lambda$ tends to 0. If the value of $s$ is finite, then it tends to

$$ s_0 = -\frac{\hbar}{\bar{\eta}} \quad (51) $$

If the value of $s$ is infinite, it is not clear which sign it has.

Next, we examine the influence of the inertia and coupling terms. More specifically, we will perform a parametric investigation of the dispersion relation between the Lyapunov coefficient, $s$, and the wavelength, $\lambda$, for different values of the inertia parameter, $\bar{\rho} = 0.001$ and $\bar{\rho} = 0.1$.

Furthermore, for each value of the inertia parameter, we investigate each coupling parameter separately, and plot the resulting dispersion relations for the real and imaginary parts of the Lyapunov coefficient. Finally, we comment on the resulting localization profile and the expected numerical response.
<table><caption>**Table 2**
Synopsis of the mathematical results concerning a multiphysics coupled rate dependent Cauchy continuum without inertia.</caption><thead><tr><th>Term</th><th>λ → 0</th><th>λ → ∞</th><th>Note</th></tr></thead><tbody><tr><td>āU</td><td>-ħ/ħ̄<br>-∞</td><td>-ħ/ħ̄<br>0</td><td>Positive maximum present, the sign of āU affects only the imaginary part.</td></tr><tr><td>āV</td><td>-ħ/ħ̄<br>-∞</td><td>Eq. (45)<br>Eq. (45)</td><td>Positive global maximum present for negative values of āV. Imaginary part not present in this case.</td></tr><tr><td>āW</td><td>-ħ/ħ̄<br>-∞/(1 + āW)</td><td>Eq. (49)<br>Eq. (49)</td><td>Positive global maximum present for negative values of āW that are larger than -1 at λ → 0.</td></tr><tr><td>A, āU</td><td>-ħ/ħ̄<br>-∞</td><td>-ħ/ħ̄<br>0</td><td>Positive local and global maximum present due to āU ≠ 0, global maximum value affected by A.</td></tr><tr><td>B, āW</td><td>-ħ/ħ̄<br>-∞/(1 + āW)</td><td>Eq. (50)<br>Eq. (49)</td><td>Strong dependence on signs. For suitable values of āW the difference lies in the value of the maximum as λ → ∞.</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** $\bar{\rho} = 0.001$. First row: Real part of the Lyapunov coefficient $Re(s)$ w.r.t. the wavelength $\lambda$. Second row: Imaginary part of the Lyapunov coefficient $Im(s)$ w.r.t. the wavelength $\lambda$. In (a,b,c,d,e) the maximum $Re(s)$ lies at $\lambda = 0$, strain localization on a mathematical plane (zero-thickness shear band) and mesh dependence develops. In (a,b,c,d), $Im(s) = 0$ at $\lambda = 0$, stationary strain localization occurs. In (e) the $maxRe(s)$ lies at $\lambda = 0$ where $Im(s) \neq 0$, traveling strain localization. In b for intermediate values of $āV$, $maxRe(s)$ lies at $\lambda \neq 0$, strain regularization occurs.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** $\bar{\rho} = 0.1$. First row: Real part of the Lyapunov coefficient $Re(s)$ w.r.t. the wavelength $\lambda$. Second row: Imaginary part of the Lyapunov coefficient $Im(s)$ w.r.t. the wavelength $\lambda$. The response stays qualitatively the same as in the case of low inertia number $\bar{\rho} = 0.001$.</figcaption></figure>

### 5.2.1. Dependence of the diffusion coefficient on the strain (āU, āV = 0, āW = 0, A = 0, B = 0)

In Figs. 4(a), 5(a), we present the real and imaginary parts of the Lyapunov coefficient $s$ versus the wavelength $\lambda$ for different values of the inertia parameter $\bar{\rho} = 0.001$, $\bar{\rho} = 0.1$.

More specifically, we see that for both values of inertia $\bar{\rho} = 0.001$, $\bar{\rho} = 0.1$ the Lyapunov coefficient consists of a real and an imaginary part. Moreover, a global maximum is present as the wavelength tends to

zero. This means that strain localization on a mathematical plane and mesh dependence of the numerical results will take place.

Furthermore, from the real part of the Lyapunov coefficient we see that next to the global maximum a local maximum is present. This has the effect that for coarse meshes, whose characteristic dimension is larger than the wavelength at the position of the local maximum, the numerical results will appear to falsely converge to a mesh independent response. However, once the characteristic mesh dimension becomes
smaller than the wavelength of the local maximum, mesh dependence of the numerical solution will ensue.

### 5.2.2. Dependence of the source term in the advection-diffusion equation on the strain ($\bar{a}V$, $\bar{a}U = 0$, $\bar{a}W = 0$, $A = 0$, $B = 0$)

In this section we study for the inertia parameters $\bar{\rho} = 0.001$, $\bar{\rho} = 0.1$ the stability of the homogeneous solution in the case where the source term of the diffusion equation depends on the strain $V$. We present in Figs. 4(b), 5(b) the values of the Lyapunov coefficient, $s$, w.r.t the wavelength, $\lambda$, for the different values of the inertia parameter.

In particular, we observe in Figs. 4(b), 5(b), that the real part of the Lyapunov coefficient presents both a real and an imaginary part. Furthermore, we note that the solution of the linearized system of (15), (17) depends on the sign and magnitude of the strain sensitivity $V$. More specifically, we note that in the case of negative coupling, $\bar{a}V < 0$, a global maximum of the Lyapunov coefficient is present for an intermediate finite value of the wavelength $\lambda > 0$. Thus, when $\bar{a}V < 0$, strain localization on an intermediate wavelength will take place. Furthermore, due to the presence of the non zero imaginary part of the Lyapunov coefficient, a traveling band of finite thickness will develop inside the layer.

Next, we focus on the positive values of the coupling parameter $\bar{a}V$. In this case we notice that the magnitude of the inertia parameter is important. More specifically, we notice that for small values of the inertia parameter $\bar{\rho} = 0.001$ (see Fig. 4(b)b) the real part of the Lyapunov coefficient presents a global maximum when the wavelength, $\lambda$, tends to zero. Strain will localize on a mathematical plane and mesh dependence will be present in the numerical results. The imaginary part of the Lyapunov coefficient is then zero, which implies that a stationary band of zero thickness will develop.

We examine next the larger values of the inertia parameter $\bar{\rho} = 0.1$, see Fig. 5(b). For positive values of the coupling parameter $\bar{a}V > 0$, we show that the Lyapunov coefficient presents a global maximum for large values of the wavelength $\lambda \to \infty$. This means that strain regularization will take place and the corresponding numerical results will be mesh independent.

### 5.2.3. Dependence of the source term in the advection-diffusion equation on the strain rate ($\bar{a}W$, $\bar{a}U = 0$, $\bar{a}V = 0$, $A = 0$, $B = 0$)

In this section we present in Figs. 4(c), 5(c) the Lyapunov coefficient in the case where the source term in (17) depends solely on the strain rate. We present in Figs. 4(c), 5(c), the analysis results for two values of the inertial parameter $\bar{\rho} = 0.001$, $\bar{\rho} = 0.1$ for varying values of the diffusion parameter $\bar{a}W$. For all values of the parameter $\bar{a}W$ and inertia $\bar{\rho}$ we examined, the values of the Lyapunov coefficient exhibit qualitatively the same characteristics.

We note that for positive values of the coupling parameter $\bar{a}W$, the real part of Lyapunov coefficient presents a finite global maximum value when the wavelength tends to zero ($\lambda \to 0$). Moreover, no imaginary part is present in this case. In the context of numerical analyses this indicates that strain will localize on the smallest mesh dimension and mesh dependence of the numerical results will take place.

Next, we investigate the Lyapunov coefficient for negative values of the coupling parameter. We note that the real part of Lyapunov coefficient tends to $+\infty$ when the wavelength becomes zero. Moreover, we note that a second local maximum is present for intermediate values of the wavelength $\lambda$. Furthermore, we note that the local maximum of the Lyapunov coefficient coincides with a region of non zero imaginary part.

In the context of numerical analyses this indicates that for coarse meshes the analysis response will exhibit a traveling strain localization profile. However, when the dimensions of the finite elements become smaller than the wavelength at the local maximum of the Lyapunov coefficient, $s$, strain will localize on the smallest mesh dimension and mesh dependence of the numerical results will occur.

### 5.2.4. Effect of advective flow and dependence of the diffusivity on the strain and $\theta$ ($A$, $\bar{a}U = 1$, $\bar{a}V = 0$, $\bar{a}W = 0$, $B = 0$)

We present in Figs. 4(d), 5(d), the analysis results for two values of the inertial parameter $\bar{\rho} = 0.001$, $\bar{\rho} = 0.1$ for varying values of the diffusion parameter $A$. For all values of the parameter $A$ and inertia $\bar{\rho}$ we examined, the values of the Lyapunov coefficient exhibit qualitatively the same characteristics. Namely, the real part of the Lyapunov coefficient tends to $+\infty$ when the wavelength tends to zero. Moreover, in the region near zero wavelength, $\lambda \to 0$, the imaginary part of the Lyapunov coefficient becomes negative. Thus, strain will localize on a traveling mathematical plane.

Moreover, we note in Figs. 4(d), 5(d), the existence of a local maximum of the real part of the Lyapunov coefficient at intermediate values of the wavelength. Its wavelength coincides with a non zero imaginary part of the Lyapunov coefficient.

In a numerical analysis we expect that for coarse meshes strain will localize on a traveling pulse of finite thickness. Moreover, for finer discretizations strain will localize on the smallest mesh dimension (strain localization on a mathematical plane).

### 5.2.5. Effect of divergence of the advective flow and dependence of the source term for $\theta$ on itself ($B$, $\bar{a}U = 0$, $\bar{a}V = 0$, $\bar{a}W = -0.5$, $A = 0$)

We present in Figs. 4(e), 5(e), the Lyapunov coefficient for two values of the inertial parameter $\bar{\rho} = 0.001$, $\bar{\rho} = 0.1$ for varying values of the advection parameter $B$. We note that for all values of the parameter $B$ and inertia $\bar{\rho}$ we examined, the values of the Lyapunov coefficient exhibit qualitatively the same characteristics. Namely, the real part of the Lyapunov coefficient tends to $+\infty$ when the wavelength tends to zero. However, in the region near zero wavelength, $\lambda \to 0$, the imaginary part of the Lyapunov coefficient does not vanish. Thus, strain will localize on a traveling mathematical plane.

Moreover, we note in Figs. 4(e), 5(e), the existence of a local maximum of the real part of the Lyapunov coefficient at intermediate values of the wavelength. Its wavelength coincides with a zero imaginary part of the Lyapunov coefficient.

In a numerical analysis, we expect that for coarse meshes strain will localize on a stationary band of finite thickness. However, for finer discretizations a traveling strain pulse will develop, which localizes on the smallest mesh dimension (traveling strain localization on a mathematical plane).

## 6. Discussion of the mathematical results

The effect of multiphysics couplings on the regularization of mechanical softening has been investigated previously by means of linear stability analysis by the authors [35] and in [42]. A direct comparison to the first work is not attempted here, as in that case simplifications in the elasto-viscoplastic response and in the dependencies of the advection-diffusion equation were assumed, which are lifted here.

[42] follow a somewhat different approach, with the coupling affecting the strength of the material only, so that the form of the determinant cannot be directly juxtaposed. A qualitative comparison of results for a similar parameter set is however possible.

More specifically, in the analysis of [42] $A = B = U = V = 0$. Moreover, the term corresponding to $W$ is set equal to one, but is multiplied by the plastic strain rate instead of the total strain rate that was considered here. For the qualitative comparison attempted here we set also $W = 1$. With the same reasoning we selected to set $\bar{a} = h_\gamma$ for the comparison, where $h_\gamma$ is the effective softening modulous defined in [42].

The results of the comparison are shown in Fig. 6. In red are plotted the results from [42], while in black those of the present analysis. It can be seen that the agreement is good, with differences being quantitative rather than qualitative. Such deviations are expected, since the formulations used for the coupling differ slightly, as discussed above.
$$ (a) \quad \bar{\alpha}=-1.0, \quad \bar{h}=0.5. $$

<figure><img src="image_9.png" /></figure>

$$ (c) \quad \bar{\alpha}=0.2, \quad \bar{h}=0.1. $$

<figure><img src="image_10.png" /></figure>

$$ (b) \quad \bar{\alpha}=-0.2, \quad \bar{h}=-0.1. $$

<figure><img src="image_11.png" /></figure>

$$ (d) \quad \bar{\alpha}=0.2, \quad \bar{h}=-0.1. $$

<figure><img src="image_12.png" /><figcaption>**Fig. 6.** Comparison between present work (in black) and [42] (in red), real parts. The two studies agree qualitatively in all cases, including positive and negative values of the coupling term, $\bar{\alpha}$, in the strain relation and strain softening $\bar{h} < 0$ vs strain hardening $\bar{h} > 0$, with the differences owing to the form of the balance equations we study.</figcaption></figure>

<caption>**Table 3**</caption>

Material parameter values for the fault geomaterial under shear. Both a dilating $\tan \psi > 0$ and a contracting material $\tan \psi < 0$ have been taken into account in order to control the sign of the coupling term $(\bar{\alpha}W)$. The scalar quantity $\alpha$ denotes the magnitude of the spherical part of the coupling tensor $\alpha_{ij}$.

<table><thead><tr><th>Parameters</th><th>Values</th><th>Units</th><th>Parameters</th><th>Values</th><th>Units</th></tr></thead><tbody><tr><td>K</td><td>666</td><td>MPa</td><td>tan ψ</td><td>-0.15 || 0.15</td><td>-</td></tr><tr><td>G</td><td>1000</td><td>MPa</td><td>H</td><td>-25</td><td>MPa</td></tr><tr><td>β<sup>s</sup></td><td>8.2 · 10<sup>-6</sup></td><td>MPa<sup>-1</sup></td><td>c<sub>hy</sub></td><td>1.0</td><td>mm<sup>2</sup>/s<sup>2</sup></td></tr><tr><td>c</td><td>100.0</td><td>MPa</td><td>θ<sub>0</sub></td><td>66.67</td><td>MPa</td></tr><tr><td>σ<sub>n</sub></td><td>200</td><td>MPa</td><td>α</td><td>-1</td><td>-</td></tr><tr><td>tan φ</td><td>0.3</td><td>-</td><td>āW</td><td>-0.63 || 0.63 · 10<sup>-2</sup></td><td>-</td></tr><tr><td>ηf<sub>0</sub></td><td>2.5 · 10<sup>-5</sup></td><td>s</td><td>H</td><td>1</td><td>mm</td></tr><tr><td>u<sub>0</sub></td><td>0.14</td><td>mm</td><td></td><td></td><td></td></tr></tbody></table>

## 7. Numerical example

In this section we aim to show numerically the localization regimes proposed in Section 5. In order to do this we focus our attention in an application of fault mechanics. In particular, we investigate the frictional response of a seismic fault under slow slip, taking into account the influence of the pore fluid pressure (see also [58]). Under such conditions, shear will localize in a thin region that we call a fault gouge (see [72–74]). The thickness of this region is related to the energy budget of the earthquake and in particular to the dissipated energy of the fault and the released energy to the surface. Correctly accounting for the friction and the dissipated energy during shearing of a fault is important for many applications including earthquake control [62,63].

### 7.1. Model description

We consider here the example of a fault gouge layer under shear, in the presence of hydro-mechanical couplings (see Fig. 1). The material properties of the fault for this specific application are shown in Table 3.

In this 1D application we consider a layer of fully saturated geomaterial inside the fault, subjected to a seismic slip shear displacement

under drained boundary conditions. We set the coupling tensor in (10) to the form $\alpha_{ij} = \alpha\delta_{ij}$, where $\alpha$ is a scalar quantity (see Table 3) and $\delta_{ij}$ is the kronecker delta. Thus, we obtain that the spherical part of the effective stress tensor depends on the pore fluid pressure $\theta$. Moreover, we set the magnitude of the scalar $\alpha$ equal to -1 in order to obtain the Terzaghi principle, which is suitable for the description of geomaterial inside the fault. Furthermore, we consider that the fault consists of an elasto-viscoplastic geomaterial. For a geomaterial the pressure dependent yield criterion and plastic multiplier are given by

$$ F(\sigma_{12}, \sigma_{22}, \varepsilon^q) = \sigma_{12} - \sigma_{22} \tan \phi - c - H \varepsilon^q, $$

$$ Q(\sigma_{12}, \sigma_{22}, \varepsilon^q) = \sigma_{12} - \sigma_{22} \tan \psi - c - H \varepsilon^q, $$

$$ \varepsilon^q = \sqrt{2e_{ij}^{vp} e_{ji}^{vp}}, \quad e_{ij}^{vp} = \varepsilon_{ij}^{vp} - \frac{1}{3} \delta_{ij} \varepsilon_{kk}^{vp} $$

$$ \dot{\varepsilon}_{ij}^{vp} = \dot{\lambda} \frac{\partial Q}{\partial \sigma_{ij}}, \qquad (52) $$

where $Q$ is the plastic potential, $c$ is the cohesion of the geomaterial, $\psi$ is the dilatancy angle and $\varepsilon^q$ is the accumulated viscoplastic deviatoric strain.

During slip, the fault gouge yields and due to the non zero dilatancy angle ($\tan \psi \neq 0$), the volume of the solid skeleton changes, introducing a pressure difference to the pore fluid. Thus, a two way coupling between the stress and mass balance equations of the problem is formed. In our example the coupling in the mass balance diffusion equation lies in the source term due to the change of the volumetric strain rate $(\bar{\alpha}W \neq 0)$. Introducing the normalized variables, $\hat{x} = \frac{x}{H}$, $\hat{u} = \frac{u}{u_0}$, $\hat{\theta} = \frac{\theta}{\theta_0}$, $\hat{t} = \frac{t}{t_0}$, where H, $u_0$ are the characteristic lengths of the problem, $\theta_0$ is the characteristic pore fluid pressure and $t_0 = \frac{H}{V_0}$ is the problem’s characteristic time dependent on the loading velocity $V_0$, the system of equations to be studied is given by

$$ \begin{aligned} \bar{h}\bar{u}_{1,22} + \bar{\eta}\dot{\bar{u}}_{1,22} + \bar{\alpha}\bar{\theta}_{,2} &= 0, \\ -W\dot{\bar{u}}_{1,2} + \dot{\bar{\theta}} - D\bar{\theta}_{,22} &= 0, \end{aligned} \qquad (53) $$

where for ease of notation we drop the hat symbol, $\hat{\cdot}$, from the normalized quantities.
We note that the system of (53) is a special case to eqs. ((15), (19)), and for the following choice of parameters

$$ 
\begin{aligned}
\bar{\eta} &= \eta f_0 \frac{V_0}{\Gamma H}, \quad \bar{h} = \frac{\Gamma - 1}{\Gamma} \frac{u_0}{H}, \quad \bar{\alpha} = -\frac{\alpha \tan \phi \theta_0}{2G\Gamma}, \\
A &= 0, \quad B = 0, \quad U = 0, \quad V = 0, \\
W &= -\frac{1}{\beta^*} \frac{\tan \psi}{K^*} \frac{u_0}{H\theta_0}, \quad D = \frac{c_{hy}}{K^*} \frac{1}{V_0 H},
\end{aligned}
\quad (54) $$

our numerical analysis corresponds to Section 5.1.3 of the mathematical analysis.

We note that in this special case, the fourth term in (15) is missing from the system (53). This has to do with the form of the material coupling tensor $\bar{\alpha}_{ij}$ that is required by the Terzaghi principle, see also Appendix A).

Furthermore, in the parameter definition of (54), $c_{hy}$ is the hydraulic diffusivity of the fault and the additional quantities $\Gamma, K^*, M$ result from the plasticity considerations in Appendix A,

$$ 
\begin{aligned}
\Gamma &= \frac{H}{2G} \sqrt{2 \left(1 + \frac{4}{9} \tan^2 \psi \right) + 1}, \\
K^* &= \left(1 - \frac{\alpha}{\beta^* M}\right).
\end{aligned}
\quad (55) $$

We note that in (53) we ignored the inertia of the fault gouge. This is because it has been shown (see [65,71]), that for typical parameter values of the fault gouge material and seismic slip velocities the inertia number becomes negligible. Furthermore, in this special case we consider that the coupling tensor, $\alpha_{ij}$, of (10) is rate independent with a zero deviatoric and a non zero spherical part (see also Appendix A).

In order to confirm the mathematical results of sub Section 5.1.3, we ran a series of numerical analyses using the Numerical Geolab (nGeo) software (see [75]), which allows for fast and efficient implementation and solution of the inelastic multiphysics numerical problem at hand.

Our analysis consists of a two step loading procedure. Initially we load the sample under isotropic normal total stress $\sigma_n = 200$ MPa and pore fluid pressure $\theta = 66.66$ MPa. Then simple shear occurs under constant displacement rate.

## 7.2. Numerical results

In Section 5.1.3 we have shown mathematically that in the presence of mechanical softening and a two way coupling in the force and diffusion equations, strain localization or strain regularization may take place, depending on the sign and magnitude of the coupling quantity $\bar{a}W$. We study two different cases of the coupling term corresponding to two different values of the parameter $\bar{a}W$, including: a) $-1 < \bar{a}W = -0.63 \cdot 10^{-2} < 0$ and b) $0 < \bar{a}W = 0.63 \cdot 10^{-2}$ (see Figs. 7 and 8 below).

### 7.2.1. Case $-1 < \bar{a}W = -0.63 \cdot 10^{-2} < 0$:

We focus first in the case of $-1 < \bar{a}W = -0.0063 < 0$. We present in Fig. 7 the numerical results for the localization width profiles of the plastic strain rate, $\dot{\lambda}$ and the corresponding shear stress-displacement diagrams, $\sigma_{12} - u_1^{top}$, for an increasing number of elements, $N$. We note that as we increase the number of elements, the strain localization width and the corresponding frictional weakening, i.e. the negative slope of the shear stress-displacement diagram, remain constant. The numerical results are mesh independent as we expect by our mathematical analysis of sub Section 5.1.3.

Considering the width of the strain localization, we can see based on the plastic strain multiplier Fig. 7, which is indicative of the plastic strain under monotonous shear, that slip is primarily concentrated in two shear bands near the boundaries of the gouge. A localization width of $\approx 0.15$ [mm] is inferred for each band, which is close to the typical thickness of the PSZ ($\approx 0.01 \sim 0.1$ [mm]) observed in nature (see for instance [74,76]).

Strain regularization occurs because a localized contractant band introduces a local increase in the rate of pore fluid pressure, that

results in a faster rate of frictional weakening in this region. Due to the equilibrium condition of the homogeneous profile, pressure in the sample needs to adjust such that the Terzaghi principle is also satisfied. This necessarily leads to a diffused pore fluid pressure increase in the neighboring cells, thus the frictional weakening is diffused, leading to strain regularization in the sample (see also [45]).

### 7.2.2. Case $0 < \bar{a}W = 0.63 \cdot 10^{-2}$:

Next, we study the case of $0 < \bar{a}W = 0.63 \cdot 10^{-2}$. We present in Fig. 8 the numerical results for the localization width profiles of the plastic strain rate, $\dot{\lambda}$, and the corresponding shear stress-displacement diagrams, $\sigma_{12} - u_1^{top}$, for an increasing number of elements $N$. We note that as we increase the number of elements, the strain localization width decreases and the corresponding frictional weakening, i.e the negative slope of the shear stress-displacement diagram also decreases. The numerical results are mesh dependent as we expect by our mathematical analysis of sub Section 5.1.3.

## 8. Conclusions

In this paper we investigated the regularization potential of a physics advection-diffusion process in the profiles of strain localization in a physics coupled mechanical system. We perform our mathematical analysis using the original form of the stress balance and advection-diffusion equations and studying the influence of all coupling terms in the system of equations.

More specifically, using the first method of Lyapunov and linearizing the system of partial differential equations around the initial homogeneous solution, we examined the localization profiles that may occur as the homogeneous solution becomes unstable. We focus on the cases where the width of the localized solution becomes zero (i.e. a zero-thickness shear band is formed). In this case the solution localizes on a mathematical plane and this translates to mesh dependent results in the numerical analyses.

Our mathematical analysis is done in several steps. Starting from the case of a strain softening strain-rate hardening material we derive the dispersion curves for the Lyapunov coefficient, $s$, and wavelength $\lambda$. We take into account the case of no inertia and the case of inertia. Our results show that in the presence of inertia and viscosity, strain will localize on a mathematical plane and mesh dependence will occur in the numerical analyses (see also [77]).

Next, we introduce the multiphysics process of advection-diffusion of a physical quantity. Again our analyses take into account the case of no inertia and the inertia case respectively. In the absence of inertia the degree of the characteristic polynomial of the dispersion relation is small and allows us to investigate in details the influence of each term of the physics process in the solution profile of the coupled system. In the presence of inertia, the degree of the characteristic polynomial is high and a full analysis of the dispersion relation is overwhelming. Thus our analysis of strain localization is done directly on the dispersion plots $s - \lambda$ for different values of the inertia parameter.

Our results show that the occurrence of strain regularization or strain localization on a mathematical plane (zero-thickness shear band) depends on the different values among the viscosity, inertia and multiphysics coupling parameters.

Considering the cases where strain regularization takes place, these indicate that specific combinations of viscosity and multiphysics parameters regularize strain and are responsible for the finite thickness of the strain profile (see also [42,45,72] among others).

Considering the cases where strain localization on a mathematical plane (zero-thickness shear band) takes place, we show that it depends on the strength of the mechanical coupling on the source term of the physics process. As a result, the addition of multiphysics is not always sufficient for regularizing the problem and mesh dependent numerical results can still take place. This is an important point for applications.
<figure><img src="image_13.png" /><figcaption>**Fig. 7.** Mesh independent analyses with parameters, $\beta^* = 8.2 \cdot 10^{-6}[MPa]^{-1}$, $\tan\psi = -0.15$, $-1 < \bar{a}W = -0.63 \cdot 10^{-2} < 0$. The curves correspond to number of elements: $N = 1, 10, 40, 80, 160, 1000$. (a) The warped displacement profile of the 1D layer at the end of each analysis. (b) The localization width of the plastic strain multiplier $\lambda$ remains constant as we increase the number of elements. (c) The negative slope of the stress displacement diagram remains constant as we increase the number of elements.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 8.** Mesh dependent analyses with parameters, $\beta^* = 8.2 \cdot 10^{-6}[MPa]^{-1}$, $\tan\psi = 0.15$, $\bar{a}W = 0.63 \cdot 10^{-2} > 0$. (a) The warped displacement profile of the 1D layer at the last converged increment of each analysis. (b) The localization width of the plastic strain multiplier $\lambda$ is always two element wide. (c) The negative slope of the stress displacement diagram decreases as we increase the number of elements.</figcaption></figure>

Our mathematical analysis on all terms of the diffusion equation of the physical quantity is general and encompasses previous analyses that studied special cases exhibiting strain localization under multiphysics couplings, see also [19,38–42,44,45,78]. We summarize our findings and comment on the occurrence of mesh dependence in Appendix C.

Finally, we showcase the results of our theoretical analysis with the help of a numerical example. More specifically, we considered a 1D numerical application in the case of the seismic slip of a mature fault gouge under the influence of pore fluid pressure. We show that our numerical and mathematical analyses agree on the range of values for the coupling strength (material properties) that results in strain localization on a mathematical plane.

From our mathematical analysis and numerical example it is clear that the use of viscosity together with a physics coupling in order to regularize strain localization and avoid mesh dependence in the numerical analyses is not always appropriate. In essence, strain regularization by the introduction of multiphysics couplings strongly depends on the material properties of the coupled system and does not constitute a general strain regularization method that leads to mesh independent numerical results.

### CRediT authorship contribution statement

**Eleni Gerolymatou:** Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing, Conceptualization, Supervision, Validation. **Alexandros Stathas:** Methodology, Software, Visualization, Writing – original draft, Writing – review & editing, Formal analysis, Investigation. **Ioannis Stefanou:** Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Supervision, Validation, Visualization, Writing – review & editing.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

No data was used for the research described in the article.

### Acknowledgements

I.S. and A.S. would like to acknowledge the support of the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program “Controlling earthQuakes - CoQuake” (Grant agreement no. 757848) and the support of the Pays de la Loire and of Nantes Metropole for the Connect Talent project “Controlling Extreme EVents - CEEV”.

## Appendix A. Numerical application: Linearized equations

### A.1. Elasto-viscoplastic behavior of a simple geo-material

We consider a general coupling form for the mechanical behavior of the geomaterial with the help of the coupling tensor, $\alpha_{ij}$ connecting the mechanical stresses and the physical quantity present in the diffusion equation. We note that this tensor can be expanded into a spherical and a deviatoric part, $\alpha_{ij} = \alpha^S \delta_{ij} + \alpha'_{ij}$, $\alpha^S = \frac{\alpha_{ii}}{3} = \alpha$, $\alpha'_{ij} = \alpha_{ij} - \alpha^S \delta_{ij}$. The coupling between the effective material stress $\sigma'_{ij}$ and the physical quantity $\theta$ is of the form

$$ \sigma_{ij} = \sigma'_{ij} + \alpha_{ij} \theta. \qquad (A.1) $$

In the case of an elastoplastic material behavior the effective stresses $\sigma'_{ij}$ are determined by the elastoplastic material matrix. Therefore, the total stress relation can be given by

$$ \sigma_{ij} = M_{ijkl}^{ep} \varepsilon_{kl} + \alpha_{ij} \theta, \qquad (A.2) $$

where $\varepsilon_{kl}$ are the total material strains. The incremental form of the mechanical component of the stress strain relationship can then be given for the stresses participating in the yield criterion,

$$ \dot{\sigma}_{12} = 2G(\dot{\varepsilon}_{12} - \dot{\varepsilon}_{12}^{vp}) + \alpha'_{12} \dot{\theta}, \qquad (A.3) $$
$$ \dot{\sigma}_{22} = M(\dot{\epsilon}_{22} - \dot{\epsilon}_{22}^{vp}) + \alpha^S \dot{\theta}, \quad (A.4) $$

where M is the 1D compression modulus for the normal stress and $\dot{\epsilon}_{12}^{vp} = \dot{\lambda}$, $\dot{\epsilon}_{22}^{vp} = \dot{\epsilon}_{12}^{vp} \tan \psi$. Replacing $\dot{\lambda} = \frac{F}{\eta F_0}$ and rearranging (A.3) we obtain

$$ \dot{\epsilon}_{12} - (2G)^{(-1)} \sigma'_{12} = \frac{\sigma'_{12} - \sigma'_{22} \tan \phi - c - H \epsilon^q}{\eta F_0}, \quad (A.5) $$

where $\epsilon^q = \epsilon_{12}^{vp} \sqrt{2(1 + \frac{4}{9} \tan^2 \psi)}$ and $\epsilon_{12}^{vp} = \epsilon_{12} - (2G)^{-1} \sigma'_{12}$. The material equation for $\dot{\sigma}_{12}$ becomes

$$ \dot{\sigma}_{12} = \alpha'_{12} \dot{\theta} + 2G \left( \dot{\epsilon}_{12} - \frac{1}{\Gamma} \dot{\epsilon}_{12} + \frac{\tan \phi}{2G\Gamma} (\dot{\sigma}_{22} - \alpha^S \dot{\theta}) + \frac{\eta f_0}{\Gamma} (\dot{\epsilon}_{12} - (2G)^{-1} (\dot{\sigma}_{12} - \alpha'_{12} \dot{\theta})) \right), \quad (A.6) $$

where

$$ \Gamma = h \sqrt{2 \left(1 + \frac{4}{9} \tan^2 \psi\right) + 1} \quad (A.7) $$

and $h = \frac{H}{2G}$, $f_0 = \frac{F_0}{2G}$. Considering that the perturbation operator acts the same way as the derivative operator, the perturbed form of the stress relation is given by:

$$ \bar{\sigma}_{12} = \alpha'_{12} \bar{\theta} + 2G \left( \frac{\Gamma - 1}{\Gamma} \bar{\epsilon}_{12} + \frac{\tan \phi}{2G\Gamma} (\bar{\sigma}_{22} - \alpha^S \bar{\rho} f) + \frac{\eta f_0}{\Gamma} (\bar{\epsilon}_{12} - (2G)^{-1} (\bar{\sigma}_{12} - \alpha'_{12} \bar{\theta})) \right). \quad (A.8) $$

We note that for the applications considered here the total normal stress $\sigma_{22}$ acting on the fault remains constant during the analysis therefore, $\bar{\sigma}_{22} = 0$. Replacing into the perturbed form of the equilibrium equation in (4) we obtain

$$ \frac{(\Gamma - 1)}{\Gamma} \bar{u}_{1,22} + \bar{\alpha}^S \frac{\tan \phi}{\Gamma} \bar{\theta}_{,2} + \frac{\eta f_0}{\Gamma} (\dot{\bar{u}}_{1,22} - (2G)^{-1} \rho \ddot{\bar{u}}_1 - (2G)^{-1} \alpha'_{12} \dot{\bar{\theta}}_{,2}) = \frac{\rho}{2G} \ddot{\bar{u}}_1. \quad (A.9) $$

Considering the Terzaghi principle to hold true, i.e. the geomaterial consists of sand or clay, $\sigma_{ij} = \sigma'_{ij} + \alpha_{ij} \theta$, ($\sigma_{ij} < 0$, $\sigma'_{ij} < 0$, $\theta > 0$ in compression, $\alpha_{ij} = -\delta_{ij}$). Thus, we obtain the following, $\alpha^S = -\delta_{ij}$, $\sigma'_{ij} = 0$, leading to $\bar{\alpha}^S_{ij} = \frac{\alpha^S}{2G} \delta_{ij} = -\frac{1}{2G} \delta_{ij}$. In this case the equilibrium equation will be written as

$$ \frac{(\Gamma - 1)}{\Gamma} \bar{u}_{1,22} + \bar{\alpha}^S \frac{\tan \phi}{\Gamma} \bar{\theta}_{,2} + \frac{\eta f_0}{\Gamma} (\dot{\bar{u}}_{1,22} - (2G)^{-1} \rho \ddot{\bar{u}}_1) = \frac{\rho}{2G} \ddot{\bar{u}}_1. \quad (A.10) $$

In the case of absence of inertia the perturbed equilibrium equation is given by

$$ \frac{(\Gamma - 1)}{\Gamma} \bar{u}_{1,22} + \bar{\alpha}^S \frac{\tan \phi}{\Gamma} \bar{\theta}_{,2} + \frac{\eta f_0}{\Gamma} \dot{\bar{u}}_{1,22} = 0. \quad (A.11) $$

## A.2. Coupled hydraulic diffusion equation

We consider now the perturbed coupled hydraulic diffusion equation:

$$ \dot{\bar{\theta}} = c_{hy} \bar{\theta}_{,22} - \frac{1}{\beta^*} \dot{\epsilon}_v \quad (A.12) $$

We notice here that the rate of increase of the partial differential equation lies on the rate of volumetric deformation of the layer $\dot{\epsilon}_v$. The volumetric change of the layer is due to the volumetric increase predicted by the dilatancy angle $\psi$ of the plastic potential,

$$ \dot{\epsilon}_v = \dot{\epsilon}_{11} + \dot{\epsilon}_{22} + \dot{\epsilon}_{33}. \quad (A.13) $$

Due to the 1D geometric constraint we have that $\dot{\epsilon}_v = \dot{\epsilon}_{22}^e + \dot{\epsilon}_{22}^{vp}$.

In order to calculate $\dot{\epsilon}_{22}$ we take into account (A.4) and the fact that during the shearing phase the total normal stress remains constant $\dot{\sigma}_{22} = 0$. In this case we have that $\dot{\epsilon}_{22} = -\alpha^S \frac{\dot{\theta}}{M} + \dot{\epsilon}_{22}^{vp}$, which using the normality rule to connect $\dot{\epsilon}_{22}^{vp}$, $\dot{\epsilon}_{12}^{vp}$ leads to

$$ \dot{\epsilon}_{22} = -\alpha^S \frac{\dot{\theta}}{M} + \dot{\epsilon}_{12}^{vp} \tan \psi. \quad (A.14) $$

<figure><img src="image_15.png" /><figcaption>Fig. B.9. Perturbation for complex s.</figcaption></figure>

We note that during shearing of the geomaterial the majority of the deformation rate is plastic. We assume then $\dot{\epsilon}_{12}^{vp} \approx \dot{\epsilon}_{12}$. Rearranging the terms we obtain,

$$ \dot{\bar{\theta}} - \frac{c_{hy}}{K^*} \bar{\theta}_{,22} + \frac{\tan \psi}{\beta^*} \frac{1}{K^*} \dot{\epsilon}_{12} = 0, \quad (A.15) $$

where

$$ K^* = \left( 1 - \frac{\alpha^S}{\beta^* M} \right). \quad (A.16) $$

## Appendix B. Solution symmetry

Up to this point it has been assumed that $k$ is positive, or at least non negative. This is common practice but is not a physical requirement. In the present section the relationship between the values of $s$ retrieved for positive values of $k$ and the values of $s$ retrieved for negative values of $k$ are investigated. It is noted that only terms with imaginary coefficients are affected by the change in sign of $k$, as those are the only ones coupled to odd powers of $k$.

For the positive valued $k_0$ the corresponding complex $s_0$ can be expressed as $\rho_s \exp i\theta_s$. The determinant of the system should be equal to zero for these values:

$$ \begin{aligned} \Delta(k_0, s_0) = & \bar{\eta} \bar{\rho} \rho_s^4 e^{i4\theta_s} + (D \bar{\eta} \bar{\rho} k^2 - i A \bar{\eta} \bar{\rho} k - B \bar{\eta} \bar{\rho} + \bar{\rho}) \rho_s^3 e^{i3\theta_s} \\ & + [(D \bar{\rho} + \bar{\eta} + \bar{\alpha} \bar{\eta} W) k^2 - i A \bar{\rho} k - B \bar{\rho}] \rho_s^2 e^{i2\theta_s} \\ & + [D \bar{\eta} k^4 + i \bar{\eta} (\bar{\alpha} U - A) k^3 + (\bar{\alpha} \eta V + \bar{\alpha} W - B \bar{\eta} + \bar{h}) k^2] \rho_s e^{i\theta_s} \\ & + D \bar{h} k^4 + i (\bar{\alpha} U - A \bar{h}) k^3 + (\bar{\alpha} V - B \bar{h}) k^2 = 0 \end{aligned} \quad (B.1) $$

Its conjugate should also be equal to zero:

$$ \begin{aligned} \bar{\Delta}(k_0, s_0) = & \bar{\eta} \bar{\rho} \rho_s^4 e^{-i4\theta_s} + (D \bar{\eta} \bar{\rho} k^2 - i A \bar{\eta} \bar{\rho} (-k) - B \bar{\eta} \bar{\rho} + \bar{\rho}) \rho_s^3 e^{-i3\theta_s} \\ & + [(D \bar{\rho} + \bar{\eta} + \bar{\alpha} \bar{\eta} W) k^2 - i A \bar{\rho} (-k) - B \bar{\rho}] \rho_s^2 e^{-i2\theta_s} \\ & + [D \bar{\eta} k^4 + i \bar{\eta} (\bar{\alpha} U - A)(-k)^3 + (\bar{\alpha} \eta V + \bar{\alpha} W - B \bar{\eta} + \bar{h}) k^2] \rho_s e^{-i\theta_s} \\ & + D \bar{h} k^4 + i (\bar{\alpha} U - A \bar{h})(-k)^3 + (\bar{\alpha} V - B \bar{h}) k^2 = 0 \end{aligned} \quad (B.2) $$

$$ \Delta(-k_0, \bar{s}_0) = 0 \quad (B.3) $$

This means that if $s_0$ is a root of the determinant for a positive $k_0$, then for $-k_0$ its conjugate $\bar{s}_0$ is a root. The real part of $s$ is not affected by the sign of $k$, but the imaginary part changes sign, if the sign of $k$ changes.

Considering the form of the perturbation in (21), the same conclusion can be drawn: if $(\bar{\alpha}_1, \bar{\theta})$ is an admissible perturbation, so is its conjugate.

## B.1. Imaginary part of roots- Traveling waves of strain localization

Due to the form assumed for the perturbations, (see (21)), the imaginary part in the value of $s$ indicates the presence of oscillations
<caption>**Table C.4**</caption>

Summary of the results of the mathematical analysis. In the absence of inertia, viscosity and multiphysics lead to mesh dependent results when the coupling terms $\bar{a}V$, $\bar{a}W$ become important. Because the mechanical component in this case does not include inertia it is not possible to decide on strain localization and mesh dependence when the diffusion equation is stable. These cases are marked as *undecided (Un)* in the table. When inertia is present strain regularization takes place when the coupling term $\bar{a}U$ becomes important for all inertia values and when the coupling term $\bar{a}V$ becomes important for large inertia values.

<table><thead><tr><th>Cauchy</th><th>Coupling</th><th>Inertia, ρ̄</th><th>ħ</th><th>ħ/η̄</th><th>D</th><th>āU</th><th>āV</th><th>āW</th><th>A</th><th>B</th><th>Mesh dependence</th></tr></thead><tbody><tr><td rowspan="10">Rate independent</td><td rowspan="2">No</td><td rowspan="2">Yes</td><td>&gt;0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>No</td></tr><tr><td>&lt;0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Yes</td></tr><tr><td rowspan="4">Yes</td><td rowspan="2">No</td><td>&gt;0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>&lt;-ħ</td><td>-</td><td>-</td><td>Yes</td></tr><tr><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>&gt;-ħ</td><td>-</td><td>-</td><td>Un</td></tr><tr><td rowspan="2">Yes</td><td>&lt;0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>&lt;-ħ</td><td>-</td><td>-</td><td>Un</td></tr><tr><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>&gt;-ħ</td><td>-</td><td>-</td><td>Yes</td></tr><tr><td rowspan="2">Yes</td><td>&gt;0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>No</td></tr><tr><td>&lt;0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Yes</td></tr><tr><td rowspan="10">Rate dependent</td><td rowspan="2">No</td><td rowspan="2">Yes</td><td>&gt;0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Yes</td></tr><tr><td>&lt;0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>No</td></tr><tr><td rowspan="8">Yes</td><td rowspan="4">No</td><td rowspan="4">&lt;0</td><td rowspan="4">&lt;0</td><td rowspan="4">&gt;0</td><td>*</td><td>0</td><td>0</td><td>0</td><td>0</td><td>No</td></tr><tr><td>0</td><td>*&lt;0</td><td>0</td><td>0</td><td>0</td><td>No</td></tr><tr><td>0</td><td>*&gt;0</td><td>0</td><td>0</td><td>0</td><td>Yes</td></tr><tr><td>0</td><td>0</td><td>*&lt;-1</td><td>0</td><td>0</td><td>Yes</td></tr><tr><td rowspan="4">Yes</td><td rowspan="4">0.001</td><td rowspan="4">&lt;0</td><td rowspan="4">&lt;0</td><td rowspan="4">&gt;0</td><td>0</td><td>0</td><td>-1 &lt;*&lt; 0</td><td>0</td><td>0</td><td>No</td></tr><tr><td>0</td><td>0</td><td>*&gt; 0</td><td>0</td><td>0</td><td>Yes</td></tr><tr><td>1</td><td>0</td><td>0</td><td>*</td><td>0</td><td>No</td></tr><tr><td>0</td><td>0</td><td>-0.5</td><td>0</td><td>*</td><td>No</td></tr><tr><td rowspan="10">Yes</td><td rowspan="8">Yes</td><td rowspan="8">0.1</td><td rowspan="8">&lt;0</td><td rowspan="8">&lt;0</td><td rowspan="8">&gt;0</td><td>*</td><td>0</td><td>0</td><td>0</td><td>0</td><td>Yes</td></tr><tr><td>0</td><td>*&lt;0</td><td>0</td><td>0</td><td>0</td><td>No</td></tr><tr><td>0</td><td>*&gt;0</td><td>0</td><td>0</td><td>0</td><td>Yes</td></tr><tr><td>0</td><td>0</td><td>*&lt;0</td><td>0</td><td>0</td><td>Yes</td></tr><tr><td>0</td><td>0</td><td>*&gt;0</td><td>0</td><td>0</td><td>Yes</td></tr><tr><td>1</td><td>0</td><td>0</td><td>*</td><td>0</td><td>Yes</td></tr><tr><td>0</td><td>0</td><td>-0.5</td><td>0</td><td>*</td><td>Yes</td></tr><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>Yes</td></tr><tr><td>0</td><td>*&lt;0</td><td>0</td><td>0</td><td>0</td><td>No</td></tr><tr><td>0</td><td>*&gt;0</td><td>0</td><td>0</td><td>0</td><td>No</td></tr><tr><td>0</td><td>0</td><td>*&lt;0</td><td>0</td><td>0</td><td>Yes</td></tr><tr><td>0</td><td>0</td><td>*&gt;0</td><td>0</td><td>0</td><td>Yes</td></tr><tr><td>1</td><td>0</td><td>0</td><td>*</td><td>0</td><td>Yes</td></tr><tr><td>0</td><td>0</td><td>-0.5</td><td>0</td><td>*</td><td>Yes</td></tr></tbody></table>

in time:

$$ \bar{u} = g \exp [s_r t + i(s_i t + kx)] \qquad (B.4) $$

The perturbation takes then the form of a traveling wave, as shown in Fig. B.9.

For the sake of a clearer illustration of the concept, it has been selected here to ignore the real part of $s$. A negative value for the imaginary part of $s$ leads to a wave traveling from left to right in the figure.

It should be noted that in the derivation presented here a uniform initial state has been assumed, i.e. parameters such as $h$, $f_0$ and $\eta$ were considered not to vary in space. Should the material exhibit softening, so that the value of $f_0$ for example is locally reduced, the appearance of traveling waves would most likely be suppressed due to the resulting geometric constraints.

## Appendix C. Analysis summary

See Table C.4.

## References

[1] Wang G, Guo X, Chen J, Han P, Su Q, Guo M, et al. Safety performance and failure criteria of lithium-ion batteries under mechanical abuse. Energies 2023;16(17):6346.

[2] Gao X, Zhu Q, Gu W. An anisotropic multiphysics damage model with application to annulus fibrosus. J Biomech 2017;61:88-93.

[3] Srivastava S, Garg RK, Sharma VS, Alba-Baena NG, Sachdeva A, Chand R, et al. Multi-physics continuum modelling approaches for metal powder additive manufacturing: A review. Rapid Prototyp J 2020;26(4):737-64.

[4] Pang M-C, Yang K, Brugge R, Zhang T, Liu X, Pan F, et al. Interactions are important: Linking multi-physics mechanisms to the performance and degradation of solid-state batteries. Mater Today 2021;49:145-83.

[5] Alpak FO. A cohesive-zone model for simulating hydraulic-fracture evolution within a fully coupled flow/geomechanics-simulation system. SPE J 2021;26(01):22-43.

[6] Andrade JE, Mital U. Multiscale and multiphysics modeling of soils. Geotech Fund Address World Chall 2019;141-68.

[7] Hueckel T. Multiphysics role in instabilities in geomaterials: A review. Instabil Model Geomech 2021;1.

[8] Regenauer-Lieb K, Veveakis M, Poulet T, Wellmann F, Karrech A, Liu J, et al. Multiscale coupling and multiphysics approaches in earth sciences: Applications. J Coupled Syst Multiscale Dyn 2013;1(3):281-323.

[9] Sciarra G, Viggiani G, Collin F. Hydromechanical instabilities in geomaterials: Advances in numerical modeling and experimental techniques. Front Mech Eng 2022;8:954665.
[10] Michopoulos JG, Apetre NA, Iliopoulos AP, Steuben JC. Elasto-plasticity, damage and multiphysics effects on the behavior of adhesive step lap joints. In: International design engineering technical conferences and computers and information in engineering conference, vol. 86212, American Society of Mechanical Engineers; 2022, V002T02A032.

[11] Michopoulos J, Mast P, Badالية R, Chwastyk T, Gause L, Farhat C, et al. Material softening issues in a multiphysics virtual wind tunnel environment. In: 40th AIAA aerospace sciences meeting & exhibit. 2002, p. 1095.

[12] Quarteroni A, Lassila T, Rossi S, Ruiz-Baier R. Integrated heart—Coupling multiscale and multiphysics models for the simulation of the cardiac function. Comput Methods Appl Mech Engrg 2017;314:345–407. http://dx.doi.org/10.1016/j.cma.2016.05.031, URL https://www.sciencedirect.com/science/article/pii/S0045782516304662. Special Issue on Biological Systems Dedicated to William S. Klug.

[13] Groen D, Zasada SJ, Coveney PV. Survey of multiscale and multiphysics applications and communities. Comput Sci Eng 2014;16(2):34–43. http://dx.doi.org/10.1109/MCSE.2013.47.

[14] Bayat M, Dong W, Thorborg J, To AC, Hattel JH. A review of multiscale and multi-physics simulations of metal additive manufacturing processes with focus on modeling strategies. Addit Manuf 2021;47:102278. http://dx.doi.org/10.1016/j.addma.2021.102278, URL https://www.sciencedirect.com/science/article/pii/S2214860421004383.

[15] O'Hara K. Reaction weakening and emplacement of crystalline thrusts: Diffusion control on reaction rate and strain rate. J Struct Geol 2007;29(8):1301–14.

[16] Burlini L, Bruhn D. High-strain zones: Laboratory perspectives on strain softening during Ductile deformation. Geol Soc Lond Special Publ 2005;245(1):1–24.

[17] Barrera O, Tarleton E, Tang H, Cocks A. Modelling the coupling between hydrogen diffusion and the mechanical behaviour of metals. Comput Mater Sci 2016;122:219–28.

[18] Sulem J, Vardoulakis I. Bifurcation analysis in geomechanics. CRC Press; 1995.

[19] Rudnicki JW, Rice J. Conditions for the localization of deformation in pressure-sensitive dilatant materials. J Mech Phys Solids 1975;23(6):371–94.

[20] Noll W. A mathematical theory of the mechanical behavior of continuous media. Arch Ration Mech Anal 1958;2(1):197–226. http://dx.doi.org/10.1007/BF00277929.

[21] Noll W. A new mathematical theory of simple materials. Found Mech Thermodyn 1974;48:243–92. http://dx.doi.org/10.1007/978-3-642-65817-4_15.

[22] De Borst R, Sluys LJ, Muhlhaus H-B, Pamin J. Fundamental issues in finite element analyses of localization of deformation. Eng Comput 1993;10(2):99–121.

[23] Goudarzi M, Gerya T, Dintner Yv. A comparative analysis of continuum plasticity, viscoplasticity and phase-field models for earthquake sequence modeling. Comput Mech 2023;1–19.

[24] Mühlhaus H-B, Vardoulakis I. The thickness of shear bands in granular materials. Géotechnique 1987;37(3):271–83. http://dx.doi.org/10.1680/geot.1987.37.3.271, http://www.icevirtuallibrary.com/content/article/10.1680/geot.1987.37.3.271. http://www.icevirtuallibrary.com/doi/10.1680/geot.1987.37.3.271. http://www.icevirtuallibrary.com/doi/10.1680/geot.1988.38.2.331.

[25] Bazant ZP. Can multiscale-multiphysics methods predict softening damage and structural failure? Int J Multiscale Comput Eng 2010;8(1).

[26] Gerard B, Pijaudier-Cabot G, Laborderie C. Coupled diffusion-damage modelling and the implications on failure due to strain localisation. Int J Solids Struct 1998;35(31–32):4107–20.

[27] Eringen AC. Microcontinuum field theories. Springer Verlag; 1999.

[28] Germain P. The Method of Virtual Power in Continuum Mechanics. Part 2: Microstructure. SIAM J Appl Math 1973;25(3):556–75. http://dx.doi.org/10.1137/0125053, URL http://epubs.siam.org/doi/abs/10.1137/0125053.

[29] Germain P. La méthode des puissances virtuelles en mécanique des milieux continus - 1. Théorie du second gradient. J Mec Theor Appl 1973;12(2):235–74.

[30] Vardoulakis I. Lecture notes on Cosserat continuum mechanics with application to the mechanics of granular media. In: 3 national meeting on "generalized continuum theories and applications". 2009, http://geolab.mechan.ntua.gr/teaching/lectnotes/CCM2009.pdf.

[31] Mindlin R. Micro-structure in linear elasticity. Arch Ration Mech Anal 1964;16(1). http://dx.doi.org/10.1007/BF00248490, URL http://www.springerlink.com/index/10.1007/BF00248490.

[32] Needleman A. Material rate dependence and mesh sensitivity in localization problems. Comput Methods Appl Mech Engrg 1988;67(1):69–85. http://dx.doi.org/10.1016/0045-7825(88)90069-2.

[33] Peirce D, Shih CF, Needleman A. A tangent modulus method for rate dependent solids. Comput Struct 1984;18(5):875–87.

[34] Heeres OM, Suiker AS, de Borst R. A comparison between the Perzyna viscoplastic model and the consistency viscoplastic model. Eur J Mech A Solids 2002;21(1):1–12.

[35] Stefanou I, Gerolymatou E. Strain localization in geomaterials and regularization: Rate-dependency, higher order continuum theories and multi-physics. ALERT Doctoral Sch 2019;780:47–85.

[36] Wang W, Sluys L, De Borst R. Viscoplasticity for instabilities due to strain softening and strain-rate softening. Internat J Numer Methods Engrg 1997;40(20):3839–64.

[37] Stathas A, Stefanou I. The role of viscous regularization in dynamical problems, strain localization and mesh dependency. Comput Methods Appl Mech Engrg 2022;388:114185.

[38] Benallal A, Bigoni D. Effects of temperature and thermo-mechanical couplings on material instabilities and strain localization of inelastic materials. J Mech Phys Solids 2004;52(3):725–53. http://dx.doi.org/10.1016/S0022-5096(03)00118-2, URL https://www.sciencedirect.com/science/article/pii/S0022509603001182.

[39] Benallal A, Comi C. On localization in saturated porous continua. Comptes Rendus de l'Académie des Sci - Ser IIB - Mech 2000;328(12):847–53. http://dx.doi.org/10.1016/S1620-7742(00)01277-0, URL https://www.sciencedirect.com/science/article/pii/S1620774200012770.

[40] Dudzinski D, Molinari A. Perturbation analysis of thermoviscoplastic instabilities in biaxial loading. Int J Solids Struct 1991;27(5):601–28.

[41] Benallal A, Comi C. Perturbation growth and localization in fluid-saturated inelastic porous media under quasi-static loadings. J Mech Phys Solids 2003;51(5):851–99.

[42] Jacquey A, Rattez H, Veveakis M. Strain localization regularization and patterns formation in rate-dependent plastic materials with multiphysics coupling. J Mech Phys Solids 2021;152:104422. http://dx.doi.org/10.1016/j.jmps.2021.104422, URL https://www.sciencedirect.com/science/article/pii/S0022509621000995.

[43] Levasseur S, Charlier R, Frieg B, Collin F. Hydro-mechanical modelling of the excavation damaged zone around an underground excavation at Mont Terri Rock laboratory. Int J Rock Mech Min Sci 2010;47(4):414–25. http://dx.doi.org/10.1016/j.ijrmms.2010.01.006, URL https://www.sciencedirect.com/science/article/pii/S1365160910000110.

[44] Plassart R, Fernandes R, Giraud A, Hoxha D, Laigle F. Hydromechanical modelling of an excavation in an underground research laboratory with an elastoviscoplastic behaviour law and regularization by second gradient of dilation. Int J Rock Mech Min Sci 2013;58:23–33.

[45] Redvand LA, Jostad HP, Grimstad G, Andresen L. Strain localisation in sensitive clays: Can rate dependency provide mesh independent results? Comput Geotech 2022;145:104678.

[46] Ahmed N, Hartmaier A. Mechanisms of grain boundary softening and strain-rate sensitivity in deformation of ultrafine-grained metals at high temperatures. Acta Mater 2011;59(11):4323–34.

[47] Collin F, Kotronis P, Pardoen B. Numerical modeling of multiphysics couplings and strain localisation. Wiley; 2021.

[48] Michopoulos JG, Farhat C, Fish J. Modeling and Simulation of Multiphysics Systems. J Comput Inf Sci Eng 2005;5(3):198–213. http://dx.doi.org/10.1115/1.2031269, arXiv:https://asmedigitalcollection.asme.org/computingengineering/article-pdf/5/3/198/6934469/198.1.pdf.

[49] Lyapunov AM. The general problem of the stability of motion. Internat J Control 1992;55:531–773, URL http://dx.doi.org/10.1080/00207179208934253.

[50] Ponthot JP. Unified stress update algorithms for the numerical simulation of large deformation elasto-plastic and elasto-viscoplastic processes. Int J Plast 2002;18(1):91–126. http://dx.doi.org/10.1016/S0749-6419(00)00097-8.

[51] Wang WM, Sluys LJ, De Borst R. Interaction between material length scale and imperfection size for localisation phenomena in viscoplastic media. Eur J Mech A Solids 1996;15(3):447–64.

[52] Sluys LJ, de Borst R, Mühlhaus HB. Wave propagation, localization and dispersion in a gradient-dependent medium. Int J Solids Struct 1993;30(9):1153–71. http://dx.doi.org/10.1016/0020-7683(93)90010-5.

[53] Routh EJ. A treatise on the stability of a given state of motion: Particularly steady motion. Being the essay to which the adams prize was adjudged in 1877, in the university of cambridge. Macmillan and Company; 1877.

[54] Hurwitz A. Ueber die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt. Math Ann 1895;46(2):273–84.

[55] Nogami M, Tomozawa M. Effect of stress on water diffusion in silica glass. J Am Ceram Soc 1984;67(2):151–4.

[56] Du X, Jin L, Zhang R. Chloride diffusivity in saturated cement paste subjected to external mechanical loadings. Ocean Eng 2015;95:1–10.

[57] Aifantis EC. On the role of gradients in the localization of deformation and fracture. Internat J Engrg Sci 1992;30(10):1279–99.

[58] Rice JR. On the stability of dilatant hardening for saturated rock masses. J Geophys Res 1975;80(11):1531–6.

[59] Chen B, Barboza BR, Sun Y, Bai J, Thomas HR, Dutko M, et al. A review of hydraulic fracturing simulation. Arch Comput Methods Eng 2021;1–58.

[60] Li Y, Zhang X, Liu X, Wang L, Zhang J. Comparative analysis of the drainage effect based on improved hydraulic flushing technology. ACS Omega 2022;7(48):44310–21.

[61] Liu T, Lin B, Fu X, Zhao Y, Gao Y, Yang W. Modeling coupled gas flow and geomechanics process in stimulated coal seam by hydraulic flushing. Int J Rock Mech Min Sci 2021;142:104769.

[62] Stefanou I, Tzortzopoulos G, Gutierrez-Oribio D. Preventing earthquake instabilities and inducing controlled, slow-slip by active fluid pressure control in the vicinity of a single seismic fault. Technical report, Copernicus Meetings; 2023.

[63] Gutierrez-Oribio D, Tzortzopoulos G, Stefanou I, Pleslan F. Earthquake control: An emerging application for robust control. theory and experimental tests. IEEE Trans Control Syst Technol 2023.
[64] Vosteen H-D, Schellschmidt R. Influence of temperature on thermal conductivity, thermal capacity and thermal diffusivity for different types of rock. Phys Chem Earth Parts a/b/c 2003;28(9-11):499-509.

[65] Rice JR, Rudnicki JW, Platt JD. Stability and localization of rapid shear in fluid-saturated fault gouge: 1. Linearized stability analysis. J Geophys Res: Solid Earth 2014;119(5):4311-33. http://dx.doi.org/10.1002/2013JB010710.

[66] Badt NZ, Huber C, Hirth G, Tullis TE. The pressure-and temperature-dependence of thermal pressurization in localized faults. J Geophys Res: Solid Earth 2023;128(8). e2023JB026558.

[67] Lasaga AC. Kinetic theory in the earth sciences. Princeton, 2014, http://dx.doi.org/10.1515/9781400864874.

[68] Sulem J, Stefanou I, Veveakis E. Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure. Granul Matter 2011;13(3):261-8. http://dx.doi.org/10.1007/s10035-010-0244-1, URL http://www.springerlink.com/index/10.1007/s10035-010-0244-1.

[69] Dal M, Fabbro R. An overview of the state of art in laser welding simulation. Opt Laser Technol 2016;78:2-14.

[70] Uday M, Ahmad Fauzi M, Zuhailawati H, Ismail A. Advances in friction welding process: A review. Sci Technol Weld Join 2010;15(7):534-58.

[71] Stathas A, Stefanou I. Fault friction under thermal pressurization during large seismic-slip: Numerical analyses and extension of the model of frictional slip. Int J Mech Sci 2023;108184.

[72] Kanamori H, Brodsky EE. The physics of earthquakes. Rep Progr Phys 2004;67(8):1429-96. http://dx.doi.org/10.1088/0034-4885/67/8/R03.

[73] Lawson A, Reid H. The mechanics of the earthquake. Calif Earthq April 1906;18:451.

[74] Sibson RH. Thickness of the seismic slip zone. Bull Seismol Soc Am 2003;93(3):1169-78. http://dx.doi.org/10.1785/0120020061.

[75] Stathas A, Stefanou I. Numerical geolab: nGeo. Zenodo; 2024, http://dx.doi.org/10.5281/zenodo.10507779.

[76] Rattez H, Stefanou I, Sulem J. The importance of Thermo-Hydro-Mechanical couplings and microstructure to strain localization in 3D continua with application to seismic faults. Part I: Theory and linear stability analysis. J Mech Phys Solids 2018;115:54-76. http://dx.doi.org/10.1016/j.jmps.2018.03.004, URL http://linkinghub.elsevier.com/retrieve/pii/S0022509617309626.

[77] Stathas A, Stefanou I. The role of viscous regularization in dynamical problems, strain localization and mesh dependency. Comput Methods Appl Mech Engrg 2022;388:114185.

[78] Loret B, Prevost J. Dynamic strain localization in fluid-saturated porous media. J Eng Mech 1991;117(4):907-22.