# Hydrodynamics of Fault Gouges From Constitutive Modelling to the Physics of Friction

Filippo Masi, Itai Einav

► To cite this version:

Filippo Masi, Itai Einav. Hydrodynamics of Fault Gouges From Constitutive Modelling to the Physics of Friction. Journal of Geophysical Research: Solid Earth, 2025, 130 (e2024JB030822), pp.1-34. 10.1029/2024JB030822. [hal-05168809](https://hal.science/hal-05168809v1)

HAL Id: hal-05168809

https://hal.science/hal-05168809v1

Submitted on 18 Jul 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
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2024JB030822

**Key Points:**

* A new physics-based pathway to study the frictional rate- and state-dependent behavior of fault gouges

* The model predicts key experimental trends of the dynamics and the energy budget of clay-rich shear zones

* The model reveals the dependence of friction on shear zone thickness and normal stress, and its independence on solid fraction

**Correspondence to:**

I. Einav,
itai.einav@sydney.edu.au

**Citation:**

Masi, F., & Einav, I. (2025). Hydrodynamics of fault gouges from constitutive modelling to the physics of friction. *Journal of Geophysical Research: Solid Earth, 130*, e2024JB030822. https://doi.org/10.1029/2024JB030822

Received 25 NOV 2024
Accepted 27 JUN 2025

**Author Contributions:**

**Conceptualization:** Filippo Masi,
Itai Einav

**Data curation:** Filippo Masi

**Formal analysis:** Filippo Masi

**Funding acquisition:** Itai Einav

**Investigation:** Filippo Masi, Itai Einav

**Methodology:** Filippo Masi, Itai Einav

**Project administration:** Itai Einav

**Resources:** Itai Einav

**Software:** Filippo Masi

**Supervision:** Itai Einav

**Validation:** Filippo Masi

**Visualization:** Filippo Masi

**Writing – original draft:** Filippo Masi

**Writing – review & editing:** Itai Einav

## Hydrodynamics of Fault Gouges From Constitutive Modelling to the Physics of Friction

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

Filippo Masi¹,² ID and Itai Einav¹ ID

¹Sydney Centre in Geomechanics and Mining Materials, School of Civil Engineering, The University of Sydney, Sydney, NSW, Australia, ²Université Grenoble Alpes, Inria, CNRS, Grenoble INP, Institute of Engineering, LJK, Grenoble, France

**Abstract** The development of rate- and state-dependent friction laws offered important insights into the key physical mechanisms of the frictional behavior of fault gouges and their seismic cycle. However, past approaches were specifically tailored to address the problem of fault shearing, leaving questions about their ability to comprehensively represent the gouge material under general loading conditions. This work establishes an alternative approach for developing a physical friction law for fault gouges that is grounded on the rigour of the hydrodynamic procedure with two-scale temperatures through *Terracotta*, a thoroughly robust constitutive model for clay in triaxial loading conditions. By specifying the model for direct shearing, the approach yields an alternative friction law that readily captures the frictional dynamics of fault gouges, including explicit dependencies on gouge layer thickness, normal stress, and solid fraction. Validated against available laboratory experiments, the friction law retains the original predictive capabilities of Terracotta in triaxial conditions and explains the rate-and-state, dilatational behavior of fault gouges in direct shear conditions. Finally, when the Terracotta friction law is connected to a spring-dashpot representation of the host rock, the combined model predicts an elastic buildup precursor to the onset of and subsequent seismicity, with results closely reflecting experimental evidence and field observations. While this study focuses on clay-rich gouges, the approach and findings are expected to offer much wider implications to a variety of materials.

**Plain Language Summary** Understanding how faults in the Earth's crust move and generate earthquakes is essential for assessing seismic hazards. Friction laws in geophysics describe the resistance that fault zones effectively transmit through gouge material. While existing friction laws have provided valuable insights, these are typically engineered directly to the specific conditions of sheared faults. This study introduces a new friction law for sheared fault gouges, which is derived from first principles using an already existing description of clay materials. Originally designed to describe the complex behavior of clay under conditions different than sheared fault gouges, the model is first generalized and is then specified to study the shear dynamics of earthquake fault zones. Thanks to the new approach, the model readily considers key factors such as gouge layer thickness, stress, and solid fraction. The new friction law accurately reproduces experimental and field observations, including critical phenomena such as fault build-up prior to a first earthquake and the post seismic cycle, all while revealing the evolution of dilatancy (volume changes). Although this research focuses on clay-rich gouges, the approach can be extended to other materials, offering broader implications for understanding earthquake mechanics and fault zone processes.

## 1. Introduction

Over recent decades, experimental and modelling efforts have focused on representing the frictional behavior of faults to deepen our understanding of seismogenesis and the seismic cycle. Laboratory tests on rock, clay, and granular fault gouges have motivated the development of rate-and-state friction laws (Dieterich, 1979; Marone, 1998; Ruina, 1983; Saffer & Marone, 2003), which aim to describe key features of earthquake nucleation and occurrence. Originally developed to phenomenologically capture the frictional and healing behavior observed in laboratory tests (Dieterich, 1979), these friction laws have been instrumental in linking experimental data with the mechanics of earthquakes and faulting (e.g., Scholz, 2019). Rate-and-state frictional parameters are generally obtained through fitting the friction law to curves obtained from velocity stepping or slide-hold-slide friction experiments (Marone, 1998; Saffer & Marone, 2003). The corresponding rate-and-state friction laws are conventionally adopted to model the mechanics of earthquakes and faulting, offering a robust explanation for numerous laboratory observations (Y. Gu & Wong, 1991; J.-C. Gu et al., 1984). Combined with continuum elasticity of the surrounding rock mass, these laws have also been applied to simulate natural and induced fault

© 2025. The Author(s).
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.
behaviors such as earthquake nucleation (Lapusta et al., 2000; Scholz, 2019), aftershock sequences (T. Chen & Lapusta, 2009; Marone et al., 1995; Segall & Rice, 1995), postseismic or interseismic fault creep and healing (Hetland et al., 2010; Marone, 1998).

Nevertheless, several limitations remain. The main challenge in using rate-and-state friction laws is the empirical nature of their parameters, which do not directly correspond to specific physical mechanisms. Consequently, existing rate-and-state friction laws cannot readily explain fault responses under diversely different laboratory and natural conditions without constant adaptation of the material constants (see J. Chen & Spiers, 2016). In particular, earthquake nucleation and rupture propagation (Rice & Ruina, 1983; Tse & Rice, 1986) depend considerably on factors such as temperature, stress level and slip velocity, which complicates the calibration of rate-and-state friction laws (Bedford et al., 2022; Carpenter et al., 2014; J. Chen & Spiers, 2016). These limitations of rate-and-state friction laws stem from two main factors: (a) their empirical foundation, which is based on a rather retrospective and idealized view of contact asperity behavior, and (b) the fact they ignore the gouge volume within the fault, which does not only sustain shear sliding but also consolidation and shear-induced dilation processes (for further details, see J. Chen & Spiers, 2016; Segall & Rice, 1995).

To address these limitations, recent effort has focused on the development of models inspired by physics at the microscale that capture steady-state and transient fault friction while grounding parameter constraints on micromechanics (among others, see Barbot, 2019, 2022; J. Chen & Niemeijer, 2017; J. Chen & Spiers, 2016; Daub & Carlson, 2008; Ikari et al., 2016; Niemeijer & Spiers, 2007; Noda, 2008; Perfettini & Molinari, 2017; Pranger et al., 2022; Van den Ende et al., 2018). For instance, Daub and Carlson (2008) and Daub et al. (2010) developed a rate-and-state constitutive model rooted on the shear transformation zone theory (see Bouchbinder & Langer, 2009; Falk & Langer, 1998; Lieou & Langer, 2012) providing an effective description of dilation and compaction phenomena in the frictional behaviour of gouges. A pivotal development along this direction is the model by Chen, Niemeijer, and Spiers, also known as the CNS model, which incorporates deformation mechanisms at the microstructural level that differentiate between a localized shear band and a bulk gouge layer. With some adaptation, the model could further include normal stress effects on the frictional behavior (Yu et al., 2024). To achieve these feats, unlike phenomenological rate-and-state friction laws, the model does not ignore the dilatational shear of faults (Makedonska et al., 2011).

While micromechanically-inspired fault models such as the CNS model have made significant strides, their ability to represent the response of the same gouge material under loading conditions that differ from direct fault shear has not been tested, leaving questions about their ability to robustly represent the material. For example, a general physical model should be capable of representing the response of the same gouge material under isotropic loading, triaxial shear conditions, or even the dependence of fault friction on the reversal of the slip velocity that may happen during an earthquake. While specifically tailored for fault shearing, such cyclical conditions should engage gouge elasticity that is not accounted for in both rate-and-state and micromechanically inspired friction laws. Similarly, most existing fault friction laws tend to focus exclusively on inelastic slip behavior, ignoring the initial elastic response of the material prior to the establishment of mature shear response. Therefore, the scope of these previous models cannot be easily extended to explain questions regarding the genesis of frictional faults. Finally, it is also important to note that these models have not fully addressed how the frictional behavior depends on the thickness of the slip layer, despite experimental and in situ evidence for this effect (e.g., Bedford & Faulkner, 2021; Byerlee & Summers, 1976; Evans, 1990; Ikari et al., 2011; Lyu et al., 2019; Makedonska et al., 2011; Marone et al., 1990; Scott et al., 1994; Scuderi et al., 2014, among others).

The objective of this work is the development of a new physics-based friction law for fault gouges that is grounded upon hydrodynamic principles. The key idea of this work lies in the philosophy of the development (see Figure 1). Instead of constructing a law directly by inspecting the specific loading conditions the gouge material sustains while the fault slides, we use an existing constitutive model for clay (Wiebicke & Einav, 2024), originally developed for triaxial loading conditions that are different to those during fault sliding. This approach ensures that our results are objective and free from potential bias. The model named *Terracotta*, leverages the rigorous hydrodynamic procedure by Landau and Lifshitz (2013), which was developed originally for fluid continua. Later, this framework was extended to describe liquid crystals (De Gennes & Prost, 1993), granular media (Jiang & Liu, 2009), and partially saturated porous media (Einav & Liu, 2023), amongst other materials. Having been grounded upon hydrodynamic principles, Terracotta satisfies the first and second law of thermodynamics, momentum and mass balance laws (Einav & Liu, 2018), material objectivity, and the thermodynamic reciprocity of
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Hydrodynamics of fault gouges: from constitutive to friction law. A fault gouge (a) subjected to tectonic movement with slip velocity $v_∞$ is represented as (b) a slip layer of thickness $h$, subjected to an overburden stress $σ_n$, experiencing slip and normal velocities (known as rates), $v$ and $v_n$, that control the development of the shear stress $τ$. The internal energy of the gouge is either held in the system through elastic stress, $σ_{ij}^e$, or dissipated as a two-stage irreversible mechanism through the generation of thermal and meso-related entropy productions ($R$ and $R_m$) as functions of the plastic strain rate, $ε_{ij}^p$, and the viscous stress, $σ_{ij}^d$, both dependent on the meso-temperature, $T^m$, while $η$ is a hydrodynamic coefficient that controls the energy decay between the scales. This hydrodynamic model is attached to a spring-dashpot system (c), with a spring of stiffness $k$ and a dashpot of damping ratio $ζ$ representing the host rock, to predict the dynamics of earthquakes, in clay-rich fault gouges.</figcaption></figure>

irreversible thermodynamic fluxes (Onsager, 1931), as well as considers the physical scale separation of temperatures in heterogeneous media and the energy flow in between (Jiang & Liu, 2009).

Building on the original Terracotta model for triaxial loading conditions, we first construct the tensorial generalization of the model, while preserving all of its original physics and predictive triaxial features. Equipped with the tensorially generalized model, we move further by specifying the unique boundary conditions of an isolated fault gouge under direct shear loading. In so doing, we derive an alternative friction law that automatically preserves the model's capacity to represent volumetric material response, including during consolidation and dilatational shear. The resulting physics-based friction law encapsulates the essential features of both rate-and-state and micromechanically inspired friction laws while naturally introducing dependencies on key meaningful parameters and state variables.

The resulting *Terracotta friction law* incorporates several physical phenomena previously ignored by fault friction laws. For example, the model accurately represents the pressure- and density-dependent elasticity of gouge media, allowing it to simulate the buildup of an elastic stress, which fully lets it account for the effects of the total normal stress. The Terracotta friction law also builds on the concept of critical state in soil mechanics (Wood, 1990) and the thermodynamic reciprocity of plastic volumetric and deviatoric plastic straining (Wiebicke & Einav, 2024). The latter feature enhances the capability of the model to predict post-peak frictional behavior, providing a physically grounded explanation for dilatational slip hardening and weakening. An important feature of the friction law is the consideration of a finite fault gouge thickness as a physical quantity. Its temporal evolution directly depends on the prescribed normal stress under direct shear boundary conditions, which also influences the elastic development of frictional instabilities. Consistent with the original Terracotta model (Wiebicke & Einav, 2024), the derived friction law accounts for the influence and evolution of the porosity (or alternatively, the solid fraction) of the gouge material. Another key ingredient is the inclusion of the meso-related temperature (from here, to be called *meso-temperature*) as a non-equilibrium variable for clay, rooted on prior developments for granular media (Alaei et al., 2021; Jiang & Liu, 2009). Unlike the conventional thermal temperature, which reflects atomic-scale thermal fluctuations, the meso-temperature captures the kinetic energy of the velocity fluctuations of the clay meso-structures, including the motion of platelets and aggregates that are overwhelmingly larger than atoms, thus attributed to have their own temperature. When stressed or deformed, these mesostructures exhibit fluctuations akin to granular particles, with pores expanding, collapsing, and clusters fragmenting. This motion at the mesoscale influences the material's elasticity and viscosity, as depicted in Figure 1. Despite the advanced hydrodynamic
procedure, the Terracotta friction law reveals a simple mathematical structure (as shall eventually be captured mathematically in Equation 24).

Accounting for these state variables and physical mechanisms enables the study of earthquake nucleation and occurrence by integrating the Terracotta friction law with conventional spring-dashpot model of the host rock, see Figure 1. The resulting Terracotta-spring model naturally incorporates the evolution and dependency of the slip zone thickness and normal stress. This model is then used to simulate the full seismogenesis cycle, including the buildup of shear stress and deformations, followed by the creation of frictional instabilities and their evolution towards limit cycles. Through a detailed parametric analysis of the key physical parameters of the Terracotta friction law, we study and quantify the influence of the: (a) slip zone thickness, (b) solid fraction, (c) normal stress, and (d) far-field velocity. The results demonstrate the predictive capabilities of this friction law in accurately capturing dependencies of the frictional behavior, in agreement with natural and induced seismicity data, previous experimental studies on clay-rich fault gouges, and discrete element method numerical simulations.

The paper is organized in a way that aims to progressively unfold the physical origin and rationale underlying each quantity and parameter in Figure 1, followed by a rigorous assessment of their validity. Section 2 introduces the physical principles guiding the hydrodynamic procedure, before briefly presenting the Terracotta constitutive model, originally formulated for triaxial testing conditions. The Terracotta model is then generalized tensorially to enable full description of any conceivable loading conditions, before being reduced to the specific loading conditions of direct shearing of fault gouges. In Section 3, we demonstrate the predictive capabilities of the resulting Terracotta friction law in predicting frictional stress-strain behaviors, including the elastic buildup of the shear stress, response under velocity stepping, and performance under different normal stresses, all as have been observed in laboratory experiments of clay-rich fault gouges. After the validation of the Terracotta friction law against controlled velocity experiments, Section 4 integrates the Terracotta friction law with a spring-dashpot analogue of the host rock to investigate fault genesis and quantify the impact of the key physical parameters on earthquake dynamics. In conclusion, we discuss the broader implications of the new model, considering its hydrodynamic foundations and potential to explore non-local features for the spatial and temporal modelling of faulting.

## 2. The Model

In striving to capture the rich, rate-dependent mechanical processes in clays, an abundance of constitutive models has been developed in geomechanics, mostly founded purely on curve fitting experimental data through a variety of empirical mathematical formulations. Recently, an alternative physics-based constitutive model, called *Terracotta* (see Wiebicke & Einav, 2024), was proposed through the rigorous hydrodynamic procedure equipped with two physical temperatures separating the scale of atoms from that of the meso-structures of clay particles and aggregates. Unlike earlier constitutive models, Terracotta captures major features of the mechanical behavior of clays under triaxial loading conditions, while avoiding to rely on abstract concepts from the mathematical theory of plasticity (cf. Alaei et al., 2021; Riley et al., 2023b). The capacity of the model to describe an astounding range of phenomena is schematically represented in Figure 2, showcasing its application to both undrained triaxial and isotropic compression, which is consistent with empirical observations of critical state in soil mechanics (see insets a and b). Additionally, Terracotta accurately predicts an impressive range of rate-dependent experimental observations, such as velocity stepping effects, relaxation, and creep phenomena (cf. inset c).

Relative to past inelastic constitutive models for clay, and despite being strongly rooted in physics, Terracotta possesses a highly simple mathematical structure that motivates studies in geophysics. Here, we therefore seek to apply this model to direct fault shearing conditions. To this end, we present hereinafter the theoretical background behind the core of the hydrodynamic procedure and outline the Terracotta model, in its original form for triaxial testing conditions. We then tensorially generalize the model to accommodate any boundary and loading conditions. We then apply the generalized model to the specific case of direct fault shearing. The resulting equations establish a new friction law for clay-rich fault gouges, which we call *Terracotta friction law* in the following, as conceptualized in Figure 3.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Schematic outline of the range of behaviors captured by the Terracotta constitutive model in its original triaxial form (re-adapted from Wiebicke & Einav, 2024). The panels showcase the model's predictions for: (a) undrained triaxial compression, (b) isotropic compression, and (c) rate-dependent processes in uniaxial compression (top) and undrained triaxial compression (bottom).</figcaption></figure>

## 2.1. Background Setting

### 2.1.1. Hydrodynamics and Two-Stage Irreversibility

In the following we summarize the main hydrodynamic principles articulated by Landau and Lifshitz (2013), Jiang and Liu (2009) and Einav and Liu (2023), on which Terracotta was founded upon. As such, unlike past constitutive models for clay, the development of Terracotta (Wiebicke & Einav, 2024) has benefitted from:

1. The conservation of energy, momentum and mass;

2. Non-negativity of the entropy production;

3. Scale separation of temperatures in heterogeneous media;

4. Understanding of energy flow across the corresponding scales;

5. Thermodynamic reciprocity (Onsager, 1931) of volumetric and deviatoric plastic fluxes;

6. Stationarity of irreversible fluxes close to thermodynamic equilibria; and

7. The relationship between thermodynamic pressure and density.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Summary of the mathematical development from (a) the original Terracotta model for triaxial conditions to its (b) tensorial generalization to its (c) specification for direct fault shearing conditions, which finally reduces to (d) the resulted Terracotta friction law. The total stress in Terracotta is equal to the sum of elastic, thermodynamic, and viscous contributions. The elastic stress (in invariant form $p^e$ and $q^e$), is expressed by Equation 10a; the thermodynamic pressure ($p^T$) is defined according to Equation 2; while the viscous stress (in invariant form $p^d$ and $q^d$) is expressed by Equation 11; and finally, the plastic transport coefficients are given by Equation 12. The resulting friction law for considers the slip velocity $v$, gouge thickness $h$, meso-temperature $T_m$, and elastic strains $\gamma^e$ and $\epsilon^e$.</figcaption></figure>
A key feature of the model is the incorporation of the meso-temperature as a non-equilibrium state variable in clay. Unlike the conventional thermal temperature, which represents the kinetic energy of the fluctuating motions of (microscopic) atoms, the meso-temperature accounts for the kinetic energy of the fluctuating motions of (mesoscopic) clay aggregates and platelets that are overwhelmingly larger than the atoms. Through the two-stage irreversibility principle of Jiang and Liu (2009), the two temperatures are used to describe the energy flow within the material from the continuum scale down to the atomistic microscale, through the mesoscale of clay aggregates and platelets.

The concept of meso-temperature, and its thermodynamic conjugated meso-entropy, is designed to capture a broad spectrum of mesoscopic states, including non-affine rearrangements and dynamic responses (Alaei et al., 2021; Goldhirsch, 2008; Jiang & Liu, 2017; Liu, 2021; Maranic et al., 2021). This is in contrast with the concept of effective temperature (Bouchbinder & Langer, 2009; Falk & Langer, 1998)—derived from Edward's configurational entropy and compactivity (Edwards & Oakeshott, 1989)—which addresses only a limited subset of configurational states, typically those associated with near-equilibrium or quasi-static conditions. Accounting for the meso-temperature is particularly important for systems under transient loadings, such as those occurring within fault gouges during transitions from stick to slip motions and vice versa.

Although not as directly measurable as its thermal counterpart, the meso-temperature can be inferred by (a) analyzing the fluctuating trajectories of meso-structures (e.g., particles in sand, and platelets or aggregates in clay, cf. Maranic et al., 2021), with advanced microscopic techniques or high-resolution simulations, or by (b) relating it to kinetic pressure (e.g., Alaei et al., 2021; Goldhirsch, 2008).

Following Wiebicke and Einav (2024), the state of clay materials is characterized by bulk density, $\rho$, elastic strain, $\epsilon_{ij}^e$, and thermal and meso-related entropies, $s$ and $s_m$, which are the thermodynamic conjugates of the thermal and meso-temperatures, respectively. Accordingly, the general form of the internal energy density is defined by

$$u \equiv u(\rho, \epsilon_{ij}^e, s, s_m). \qquad (1)$$

The complete derivation of the hydrodynamic procedure can be found in Alaei et al. (2021) for materials that involve exactly the same set of state variables. Following the same procedure, the total stress is defined as

$$\sigma_{ij} \equiv \sigma_{ij}^e + \sigma_{ij}^d + p^T \delta_{ij}, \qquad \sigma_{ij}^e \equiv \frac{\partial u}{\partial \epsilon_{ij}^e}, \quad p^T \equiv -\left. \frac{\partial(u/\rho)}{\partial(1/\rho)} \right|_{\epsilon_{ij}^e, \rho^e, \rho^e} \qquad (2)$$

where $\sigma_{ij}^e$ is the elastic stress tensor, $p^T$ is the thermodynamic pressure, and $\sigma_{ij}^d$ is the dissipative viscous stress, with $\delta_{ij}$ being the Kronecker delta.

Next, we determine the evolution equations for the four state variables embodied in Equation 1. Noting that throughout the article, Einstein's summation is implied for repeated indices $i, j = 1, 2, 3$, the evolution of the bulk density is derived from the mass balance, $\dot{\rho} + \nabla_i (\rho v_i) = 0$, with $v_i$ being the velocity vector. We adopt the usual sign convention in soil mechanics that consider positive deformation in compression and identify the trace of the velocity gradient with the negative volumetric strain rate, $\nabla_i v_i = -\dot{\epsilon}_{kk} \equiv -\dot{\epsilon}_v$. By further considering dry fault gouges with negligible air mass relative to the solid mass, the solid fraction could be calculated using $\phi \equiv \rho/\rho_s$, where $\rho_s$ is the solid density defined as the ratio of the solid mass over the solid volume. The evolution equation for the solid fraction is then defined by

$$\dot{\phi} + v_i \nabla_i \phi = \phi \dot{\epsilon}_v, \qquad (3)$$

under the assumption of incompressible solid particles, that is, $\dot{\rho}_s \approx 0$. Modifications that account for air and water densities within the hydrodynamic procedure would go beyond the current scope of this paper, but are certainly possible and would follow something akin to the work in (Y. Chen & Einav, 2024; Y. Chen et al., 2023).

The evolution for the elastic strain is given by Einav and Liu (2018) and reads
$$ \dot{\epsilon}_{ij}^{e} + v_{k} \nabla_{k} \epsilon_{ij}^{e} + \Omega_{ik} \epsilon_{kj}^{e} - \epsilon_{ik}^{e} \Omega_{kj} = \dot{\epsilon}_{ij} - \dot{\epsilon}_{ij}^{p}, \quad (4) $$

where $\dot{\epsilon}_{ij}$ and $\Omega_{ij}$ are respectively the symmetric and skew-symmetric parts of the spatial velocity gradient, $\nabla_j v_i$, and $\dot{\epsilon}_{ij}^p$ is the inelastic strain rate. Following scale separation (Jiang & Liu, 2009), where the granular entropy is replaced by the more general meso-related one (Einav & Liu, 2023; Wiebicke & Einav, 2024), the thermal and meso-related entropies satisfy the following balance equations:

$$ \dot{s} + \nabla_i (sv_i - f_i) = \frac{R}{T} \geq 0, \quad (5a) $$

$$ \dot{s}_m + \nabla_i (s_m v_i - f_i^m) = \frac{R_m}{T_m}, \quad (5b) $$

where $T \equiv \frac{\partial u}{\partial s}$ and $T_m \equiv \frac{\partial u}{\partial s_m}$ are the thermal and the meso-temperatures, respectively; $sv_i$ and $s_m v_i$ are the thermal and meso-related convective entropy currents; $f_i$ and $f_i^m$ are the thermal and meso-related dissipative entropy currents; and $R$ and $R_m$ are the thermal and meso-related entropy productions. Note that only the former is required to be positive-semidefinite by virtue of the second law of thermodynamics.

By combining Equations 3–5 with the momentum and energy balances, and the principle of two-stage irreversibility in the seminal work by Jiang and Liu (2009), the total rate of entropy production comprises two-scale contributions:

$$ R + R_m \geq 0 \quad (6a) $$

$$ R_m = \sigma_{ij}^d \dot{\epsilon}_{ij} + f_i^m \nabla_i T_m - \bar{\eta} T_m^2, \quad (6b) $$

$$ R = \dot{\epsilon}_{ij}^p \sigma_{ij}^e + f_i \nabla_i T + \bar{\eta} T_m^2 \geq 0, \quad (6c) $$

where the viscous dissipation $\sigma_{ij}^d \dot{\epsilon}_{ij}$ is entirely attributed to the meso-related entropy production, unlike Jiang and Liu (2009) who placed part of it in the thermal entropy. The dissipated energy within the system can either agitate the microstructure, thereby elevating the thermal temperature, $T$, or the mesostructures, thus elevating $T_m$, while overtime the energy can flow from the meso- to the micro-scale in a two-stage process, see Figure 1b. The flow between the scales is given by the $\pm \bar{\eta} T_m^2$ sink/source terms, with $\bar{\eta} \geq 0$ so that the energy will indeed flow from the higher to the lower scale.

Next, we look at the dissipative fluxes in the entropy production Equation 6, namely: $f_i$, $f_i^m$, $\sigma_{ij}^d$, and $\dot{\epsilon}_{ij}^p$ and determine their hydrodynamic equations in terms of Onsager's reciprocity relationships (Onsager, 1931). Following (Wiebicke & Einav, 2024) these are specified by

$$ f_i = \frac{k_T}{T} \nabla_i T, \quad f_i^m = \frac{k_m}{T_m} \nabla_i T_m, \quad \sigma_{ij}^d = \mathbf{X}_{ijkl}^d \dot{\epsilon}_{kl}, \quad \dot{\epsilon}_{ij}^p = \mathbf{X}_{ijkl}^p \sigma_{kl}^e, \quad (7) $$

where $k_T$ and $k_m$ are the micro- and meso-scopic diffusivity coefficients, $\mathbf{X}_{ijkl}^p$ and $\mathbf{X}_{ijkl}^d$ are the transport coefficients, which could in general be functions of the state variables. Following Onsager (1931), the transport coefficients must fulfil a number of reciprocity requirements stemming from time-reversal symmetry and the entropy production inequality in Equation 6, as described by Wiebicke and Einav (2024).

### 2.1.2. Original Triaxial Form of the Model

Building upon the above theoretical and hydrodynamic settings, the Terracotta model (Wiebicke & Einav, 2024) was originally developed for standard triaxial conditions, in terms of triaxial stress and strain invariants. The isotropic stress and volumetric elastic strain invariants are defined as $p = \sigma_{ii}^e$ and $\epsilon_v^e = \epsilon_{ii}^p$, respectively. The
deviatoric stress and elastic strain invariants read $q = (3/2s_{ij}s_{ij})^{1/2}$ and $\epsilon_s^e = (2/3e_{ij}^e e_{ij}^e)^{1/2}$, where $s_{ij} = \sigma_{ij} - p\delta_{ij}$ and $e_{ij}^e = \epsilon_{ij}^e - \frac{\epsilon_s^e}{3}\delta_{ij}$ are the deviatoric stress and elastic strain tensors, respectively.

The original form of the model considers isothermal processes, hence the dependency of the internal energy with respect to the thermal entropy is neglected. The internal energy density, $u$, in Equation 1, is decomposed into elastic and meso-related contributions, $u^e$ and $u^m$, through:

$$u = u^e(\phi, \epsilon_{ij}^e) + u^m(s^m), \quad u^e = \phi^6 \left( \frac{\tilde{K}}{6} \epsilon_v^{e3} + \frac{3}{2} \tilde{G} \epsilon_v^{e2} \epsilon_s^e \right), \quad u^m = \frac{\Gamma}{4} s_m^2, \quad (8)$$

where $\tilde{K}$ and $\tilde{G}$ represent the intrinsic bulk and shear stiffness constants, respectively; and $\Gamma = 1 \text{ K}^2/\text{kPa}$ is a positive constant controlling the contribution of meso-related entropy to the internal energy (its inverse could be considered as a meso-related analogue of the thermal heat capacity).

By virtue of the stress decomposition in Equation 2, the volumetric and deviatoric stress invariants are given by

$$p = p^e + p^d + p^T, \quad (9a)$$

$$q = q^e + q^d, \quad (9b)$$

where $p^e, q^e$ and $p^d, q^d$ are respectively the elastic and dissipative (viscous) components. The elastic volumetric and deviatoric stresses and the thermodynamic pressure are obtained from Equation 2, giving

$$p^e \equiv \frac{\partial u^e}{\partial \epsilon_v^e} = \phi^6 \left( \frac{\tilde{K}}{2} \epsilon_v^{e2} + \frac{3\tilde{G}}{2} \epsilon_s^{e2} \right), \quad (10a)$$

$$q^e \equiv \frac{\partial u^e}{\partial \epsilon_s^e} = 3\tilde{G}\phi^6 \epsilon_v^e \epsilon_s^e, \quad (10b)$$

$$p^T = \frac{T_m^2}{\Gamma} \Leftrightarrow T_m = \frac{\Gamma}{2} s^m, \quad (10c)$$

where the elastic contribution in the thermodynamic pressure, that is, $5u^e$, was found negligible (Wiebicke & Einav, 2024) and is thus omitted for the sake of simplicity. The above equations are completed by the constitutive relationships for the dissipative fluxes, $p^d$ and $q^d$, and the evolution equations for the solid fraction, Equation 3, the meso-temperature, $T_m$, and the elastic strain, $\epsilon_v^e$ and $\epsilon_s^e$. Their derivation follows the hydrodynamic procedure, cf. Section 2.1.1, and the exact procedure is sketched below for the generalised model in tensorial form. For the original derivation in triaxial conditions, we refer to Wiebicke and Einav (2024), with the corresponding expressions presented in Figure 3a.

More precisely, the calibration of the triaxial invariants of the viscous stress, $p^d$ and $q^d$, is performed in the simplest way to ensure the semi-positivity of the viscous dissipation, leading to the following formulation:

$$p^d = aT_m\dot{\epsilon}_v, \quad q^d = cT_m\dot{\epsilon}_s, \quad (11)$$

where the viscous stress is made to decay to zero as the meso-temperature approaches thermodynamic equilibrium ($T_m \to 0$), as expected from non-equilibrium thermodynamic properties.

Next, the stationary value of the meso-temperature is used to calibrate the dependence of the plastic transport coefficients, so to enable the recovery of experimental critical state observations. This calibration approach was proposed by Wiebicke and Einav (2024) and was informed by two canonical asymptotic states: normal consolidation state emerging by continuous isotropic compression and the critical state emerging by continuous shearing (Wiebicke & Einav, 2024). Following this procedure, the following plastic transport coefficients were obtained:
$$a^p = \sqrt{\frac{\eta}{a}} \frac{1}{p_c(\phi)} \qquad b^p = -a \frac{q^e}{M^2 p^e}, \qquad c^p = \sqrt{\frac{\eta}{c}} \frac{1}{M \omega p_c(\phi)} + \frac{a}{M^2}, \qquad p_c(\phi) = p_l \phi^\lambda, \qquad (12)$$

where a, c are positive material constants, ω and M are the parameters describing the critical state line (see Figure 2 and refer to Wiebicke & Einav, 2024), and p_c(φ) describes the stress-density relationship on the isotropic compression line, where λ is the slope of the isotropic compression line on a bi-logarithmic plot of pressure against solid fraction and p_l = φ_T^λ (kPa) is calibrated with φ_l representing the solid fraction at a reference pressure of 1 kPa, see Figure 2.

## 2.2. Generalization to the Tensorial Form

Based on the hydrodynamic procedure, we extend the above form of the Terracotta model to its tensorial form with the objective of addressing more general loading conditions than just triaxial loading, including those related to the shearing of fault gouges. By following the specification of the internal energy density given in Equation 8, we generalize the elastic contribution in terms of the full elastic strain tensor, ε^e, namely

$$u^e = \phi^6 \left( \frac{\tilde{K}}{6} \varepsilon_v^{e3} + \tilde{G} \varepsilon_v^e e_{ij}^e e_{ij}^e \right), \qquad (13)$$

which is equivalent to the original expression in Equation 8 when specified in terms of triaxial invariants. As in the triaxial form of Terracotta, we neglect the dependence of the internal energy on the thermal temperature. However, such dependencies could easily be accounted for through the hydrodynamic procedure, as further discussed in paragraph Section 2.2.1.

By virtue of Equation 2, the elastic stress tensor becomes

$$\sigma_{ij}^{e} \equiv \frac{\partial u}{\partial \varepsilon_{ij}^{e}} = \phi^{6} \left( \frac{\tilde{K}}{2} \varepsilon_{v}^{e2} \delta_{ij} + \tilde{G} \delta_{ij} e_{kl}^{e} e_{kl}^{e} + 2 \tilde{G} \varepsilon_{v}^{e} e_{ij}^{e} \right), \qquad (14)$$

whose volumetric and deviatoric invariants are exactly the same as those in the triaxial formulation, that is, Equation 10. The mesoscopic component of the internal energy, u^m, is kept the same as in Equation 8, ensuring that the resulting thermodynamic pressure in Equation 10c, p^T, also remains unchanged.

It is worth noticing that the elastic energy density in Equation 13 reveals a physical elastic instability when the determinant of the instantaneous stiffness matrix, that is, the Hessian ∂σ^e/∂ε^e, vanishes. The instability criterion is equivalent with that obtained in the original triaxial form (see Wiebicke & Einav, 2024) and can be expressed in terms of the ratio between the deviatoric and volumetric elastic stresses:

$$\frac{q_e}{p_e} \leq M_e, \qquad M_e = \sqrt{\frac{3\tilde{G}}{\tilde{K}}}. \qquad (15)$$

### 2.2.1. Thermal and Meso-Related Heat Equations

The concept of two-stage irreversibility (Jiang & Liu, 2009) provides us with the form of the evolution equations for the micro- and meso-temperatures. The full hydrodynamic derivation of their corresponding rate equations can be found in Riley et al. (2023b, Appendix A), considering the two-stage irreversibility. In the absence of convective entropy currents, the evolution equations for the micro- and meso-temperatures are given by

$$\dot{T} = \frac{k_T}{c_T \rho} \nabla_i^2 T + \dot{\varepsilon}_{ij}^p \sigma_{ij}^e + \bar{\eta} T_m^2 \qquad (16a)$$

$$\frac{2T_m}{\Gamma} \dot{T}_m = k_m \nabla_i^2 T_m + \sigma_{ij}^d \dot{\varepsilon}_{ij} - \bar{\eta} T_m^2, \qquad (16b)$$
where $c_T$ is the micro-related specific thermal heat capacity. Note that the Laplace operator, $\nabla_i^2$, represents the diffusion of the thermal and meso-temperatures and can physically capture second-gradient phenomena, known as non-local effects in granular physics (Henann & Kamrin, 2013; Kamrin & Koval, 2012). Predicting these phenomena requires spatial integration of the underlying boundary value problem, which is beyond the focus of the current paper. From here on, we shall therefore neglect the second gradient terms and focus on phenomena that happen homogeneously within representative elementary volumes, independent of non-local effects.

### 2.2.2. Dissipative Fluxes and Calibration

By virtue of the time-reversal symmetry requirements, the relationships for the viscous stress and plastic strain rate introduced in Equation 7 are generalised in term of the volumetric invariant and the deviatoric tensor as follows,

$$ \begin{bmatrix} p^d \\ s^d_{ij} \end{bmatrix} = \frac{2T^m}{\Gamma} \begin{pmatrix} a^d & b^d_{kl} \\ b^d_{ij} & c^d_{ijkl} \end{pmatrix} \begin{bmatrix} \dot{\epsilon}_v \\ \dot{\epsilon}_{kl} \end{bmatrix}, \quad \begin{bmatrix} \dot{\epsilon}_v^p \\ \dot{\epsilon}_{ij}^p \end{bmatrix} = T_m \begin{pmatrix} a^p & b^p_{kl} \\ b^p_{ij} & c^p_{ijkl} \end{pmatrix} \begin{bmatrix} p^e \\ s^e_{kl} \end{bmatrix}. \quad (17) $$

Next, following Wiebicke and Einav (2024), we consider the energy sink coefficient to be linearly proportional to the meso-temperature, that is, $\bar{\eta} = 2 \frac{\eta}{T_m}$. We then substitute the constitutive relationship for the viscous stress in the meso-related heat Equation 16b, so by neglecting non-local effects we obtain the following evolution equation for the meso-temperature

$$ \dot{T}_m = a \dot{\epsilon}_v^2 + \frac{2}{3} c \dot{\epsilon}_{ij} \dot{\epsilon}_{ij} - \eta T_m^2, \quad (18) $$

where in order to retrieve the expression given in Wiebicke and Einav (2024, Equation 36), we use

$$ a^d = a, \quad b^d_{ij} = 0_{ij}, \quad c^d_{ijkl} = \frac{2}{3} c \delta_{ik} \delta_{jl}, \quad (19) $$

where $a$ and $c$ are positive material constants.

Using the above expressions and following Wiebicke and Einav's calibration procedure, we obtain the following plastic transport coefficients for the tensorially general case:

$$ a^p = \sqrt{\frac{\eta}{a}} \frac{1}{p_c(\phi)}, \quad b^p_{ij} = -\frac{3}{2} \frac{a^p}{M^2} \frac{s^e_{ij}}{p^e}, \quad c^p_{ijkl} = \frac{3}{2} \left( \sqrt{\frac{\eta}{c}} \frac{1}{M \omega p_c(\phi)} + \frac{a^p}{M^2} \right) \delta_{ik} \delta_{jl}. \quad (20) $$

Given the specific form of the transport coefficients, the full tensorial expressions for the viscous stress and the plastic strain rate become.

$$ \sigma_{ij}^d = \frac{2T^m}{\Gamma} \left( a \dot{\epsilon}_v \delta_{ij} + \frac{2}{3} c \dot{\epsilon}_{ij} \right), \quad (21a) $$

$$ \dot{\epsilon}_{ij}^p = \frac{T^m}{M^2 p_1 \phi^2} \sqrt{\frac{\eta}{a}} \left[ \left( M^2 - \left( \frac{q^e}{p^e} \right)^2 \right) \frac{p^e}{3} \delta_{ij} + \frac{3}{2} \sqrt{\frac{a}{c}} \frac{M}{\omega} s_{ij}^e \right]. \quad (21b) $$

It is worth noticing that when evaluated in terms of the volumetric and deviatoric invariants—namely, $p^d$, $q^d$ and $\dot{\epsilon}_v^p$, $\dot{\epsilon}_s^p$—the above expressions yield the exact same form of the transport coefficients as those derived in the original triaxial formulation of Terracotta model. Moreover, note that the requirement for the positive semi-definiteness of the plastic transport coefficients coincides with that identified in Wiebicke and Einav (2024, Equation 53), namely
$$ \frac{q_e}{p_e} \le M_o, \qquad M_o = M \sqrt{\frac{a}{c} \frac{M}{\omega} + 1}. \qquad (22) $$

## 2.3. Specification of the Terracotta Friction Law

From the newly expressed tensorial form of the Terracotta model, we specify the boundary conditions of an isolated fault gouge and derive the corresponding friction law under direct shearing configuration. In so doing, we use plane strain conditions and adopt the following minimal Voigt notation:

$$ \sigma_{ij} \rightarrow \begin{bmatrix} \sigma_n \\ \sigma_\ell \\ \tau \end{bmatrix}, \qquad \epsilon_{ij}^e \rightarrow \begin{bmatrix} \epsilon_n^e \\ \epsilon_\ell^e \\ \gamma^e \end{bmatrix}, $$

with superscripts n, ℓ referring to the normal and longitudinal directions, respectively. The equivalence with the tensorial form is given by $\sigma_n = \sigma_{11}$, $\sigma_\ell = \sigma_{22} = \sigma_{33}$, $\tau = \sigma_{12}$ and $\epsilon_n^e = \epsilon_{11}^e$, $\epsilon_\ell^e = \epsilon_{22}^e = \epsilon_{33}^e$, $\gamma^e = 2\epsilon_{12}^e$.

Next, we consider a fault gouge of constant length L and variable thickness h, whose state is determined by the following initial conditions at thermodynamic equilibrium

$$ h(0) = h^0, \quad \phi(0) = \phi^0, \quad T_m(0) = 0, \quad \epsilon_{ij}^e(0) = \begin{bmatrix} \epsilon_n^{e0} \\ \epsilon_\ell^{e0} \\ 0 \end{bmatrix}, \quad \sigma_{ij}(0) = \begin{bmatrix} \sigma_n^{e0} \\ \sigma_\ell^{e0} \\ 0 \end{bmatrix}, \qquad (23) $$

where the initial conditions for the elastic strain are obtained by equating the elastic stress in Equation 14 with the prescribed values for $\sigma_n^e$, $\sigma_\ell^e$, under either isotropic or uniaxial stationary compression.

Starting from the initial state in Equation 23, the fault gouge is subjected to a slip velocity v, see Figure 3d. To ensure a constant normal stress $\sigma_n = \sigma_n^0$, the material undergoes dilation or contraction along the normal direction, following the evolution of the gouge thickness h and its compaction velocity, $v_n$, again see Figure 3d. With that, the resulting Terracotta friction law reads

$$ \tau \left( \underbrace{\epsilon^e}_{\text{state}}, \gamma^e, T_m, \overbrace{h, t, v, v_n}^{\text{rate}} \right) = \frac{A}{h^6} \epsilon_v^e \gamma^e + \frac{B}{h} T_m v \qquad (24) $$

where the constants $A \equiv \tilde{G}(\phi^0 h^0)^6$ and $B \equiv (2/3)c/\Gamma$. The above equation reflects the main result of the paper, which is formulated in terms of the following evolution equations for the state variables

$$ \dot{\epsilon}_v^e = \dot{\epsilon}_n - \frac{T_m p^e}{p_c(\phi)} \sqrt{\frac{\eta}{a}} \left(1 - \left(\frac{q^e}{Mp^e}\right)^2\right), \qquad (25a) $$

$$ \dot{\gamma}^e = \dot{\gamma} - \frac{3T_m \tau^e}{M \omega p_c(\phi)} \sqrt{\frac{\eta}{c}}, \qquad (25b) $$

$$ \dot{T}_m = \left( a + \frac{4c}{9} \right) \dot{\epsilon}_n^2 + \frac{c}{3} \dot{\gamma}^2 - \eta T_m^2, \qquad (25c) $$

$$ \dot{h} = -h \dot{\epsilon}_n, \qquad (25d) $$
where $\phi = (\phi^0 h^0)/h$. The total strain rates are then given in terms of the normal and slip velocities that act as the rate variables in the present model:

$$ \dot{\gamma} = \frac{v}{h} \quad \text{such that} \quad v = v(t), \qquad (26a) $$

$$ \dot{\epsilon}_n = -\frac{v_n}{h} \quad \text{such that} \quad \dot{\sigma}_n(\text{state}, v, v_n) = 0, \qquad (26b) $$

where $v(t)$ is imposed as function of time. Appendix A presents further details of the numerical implementation of the friction law (cf. Masi & Einav, 2024).

Furthermore, Appendix B presents an analytic derivation of the steady-state friction stress $\tau_s$ and coefficient $\mu_s$, under constant slip velocity $v$ and normal stress $\sigma_n$. Unlike conventional rate-and-state friction laws that consider a purely viscous model with logarithmic dependence on slip velocity, the Terracotta friction law embodies an inherent elastic component, proportional to the critical state line, $\propto M$, and a viscous term, proportional to the square of the velocity, $\propto v^2/h_s$, where $h_s$ is the gouge thickness at steady-state. Although the quadratic nature of the steady-state friction law limits the applicability of the current simple Terracotta friction law over extremely broad velocity ranges, as opposed to empirical logarithmic-like scaling in laboratory experiments, the new law provides robust predictive capabilities within the velocity regimes considered in the present analyses. In the applications to earthquakes of the Terracotta-spring model, the steady state develops under relatively small velocities, rendering the steady state limit of the friction law to be quite realistic. Nevertheless, it is envisaged that the Terracotta friction law could be reconciled with a logarithmic-like scaling by altering the current linear scaling of the rheological coefficients (i.e., $p_c$, $a$, $b$ and $c$) on the meso-temperature $T_m$.

In the following, the shear stress evolution predicted by the Terracotta friction law is expressed in terms of the friction coefficient, that is, $\mu = \tau/\sigma_n$.

### 2.3.1. Calibration of Material Parameters

The Terracotta friction law (Equations 24–26) requires the specification of nine material parameters ($\tilde{K}$, $\tilde{G}$, $M$, $\omega$, $\lambda$, $\phi_l$, $a$, $c$, and $\eta$). However, when compared to conventional rate-and-state friction laws, which in their simplest form involve four to five parameters, the calibration of these parameters follows relatively simple and standard tests (Wiebicke & Einav, 2024) and account for much broader physics. In particular, the two elastic constants, $\tilde{K}$ and $\tilde{G}$, allow us to describe pressure and density dependent elastic responses and can be calibrated using conventional elastic wave measurements (cf. Riley et al., 2023a) or by fitting the model to observed stress-strain responses, as performed hereafter. Parameters $\phi_l$ and $\lambda$, which define inelastic compressibility that also goes beyond conventional both rate-and-state and micromechanically inspired friction laws, can be directly calibrated from isotropic/uniaxial compression tests or inferred from the post-peak slope in triaxial compression or direct shear tests. The two critical state parameters, $M$ and $\omega$, can be obtained from triaxial compression or direct shear tests, with $\omega = 0.5$ being a good value in the absence of density measurements, in agreement with the tacit assumption for clay in the modified cam clay model (Wood, 1990) of soil mechanics. The remaining three rheological parameters, $\eta$, $a$, and $c$, require calibration thorough rate-dependent experiments such as velocity stepping or slide-hold-slide tests, as those considered below.

### 2.3.2. Calibration of System Parameters

Together with the above material parameters, it is necessary to determine appropriate initial conditions for the normal stress, $\sigma_n^0$, solid fraction, $\phi^0$, and layer thickness, $h^0$. The normal stress is conventionally prescribed or measured in laboratory experiments (cf. Saffer & Marone, 2003), or in the field may be estimated as the overburden stress from the fault gouge depth (cf. Kanamori & Brodsky, 2004). The solid fraction in general varies through the thickness of the fault gouge (see Casas et al., 2023), though as discussed earlier all such gradients are neglected in the present analysis. One way to determine the initial layer thickness is to follow Rice (2006), who associated it with the localized shear zone within the gouge of less than 1 ÷ 5 mm. This shear zone of an ultracataclastic gouge is influenced by factors such as fault maturity and the chemical composition of its materials. It is worth noticing that in the present formulation, we consider a nominal uniform shear strain rate across the
<caption>**Table 1**</caption>

Parameters for the Terracotta Friction Law Calibrated for the Experiments Conducted by Saffer and Marone (2003) and Bedford et al. (2022) for Illite Shale, Smectite, and Kaolinite

<table><thead><tr><th rowspan="2">Gouge</th><th colspan="6">Mechanical parameters</th><th colspan="3">Rheological parameters</th></tr><tr><th>K̃ (MPa)</th><th>G̃ (MPa)</th><th>M (-)</th><th>ω (-)</th><th>λ (-)</th><th>φ<sub>I</sub> (-)</th><th>a (Ks)</th><th>c (Ks)</th><th>η (K<sup>-1</sup>s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Illite</td><td>1,040</td><td>400</td><td>0.96</td><td>0.6</td><td>10</td><td>0.20</td><td>7</td><td>7</td><td>7.2·10<sup>4</sup></td></tr><tr><td>Smectite</td><td>260</td><td>100</td><td>0.42</td><td>0.5</td><td>10</td><td>0.28</td><td>2.4·10<sup>4</sup></td><td>2.4·10<sup>4</sup></td><td>10</td></tr><tr><td>Kaolinite</td><td>15.6</td><td>12.5</td><td>0.44</td><td>0.4</td><td>4.5</td><td>0.05</td><td>0.42</td><td>10</td><td>3·10<sup>6</sup></td></tr></tbody></table>

gouge thickness, which neglects the coseismic evolution of shear localization from a full gouge thickness to a localized band. Nevertheless, despite precluding this important feature (cf. Barras & Brantut, 2025; Platt et al., 2014; Rice, 2006), this simplification readily enables to account for the dependence of the frictional response of the shear zone thickness, cf. Sections 3 and 4, in a straightforward way. Further accounting for the dynamic evolution of the shear zone profile by solving the underlying boundary value problem can, in future, provide deeper insights into the transition from distributed to localized slip.

In clay, on which the current paper focuses, the thickness of the localized shear zones in the fault gouge is considered to be approximately 10 to 100 times the clay platelet diameter, which is of the order of 1 µm (Bigaroni et al., 2023; Haines et al., 2013; Morgenstern & Tchalenko, 1967; Rice, 2006; Volpe et al., 2024). Although studies have explored the effects of earthquake nucleation and seismic events on the layer thickness using both laboratory samples and natural fault gouges (e.g., Bedford & Faulkner, 2021; Byerlee & Summers, 1976; Evans, 1990; Ikari et al., 2011; Lyu et al., 2019; Marone et al., 1990; Scott et al., 1994; Scuderi et al., 2014, among others), the quantification of this relationship remains elusive. Therefore, the impact of layer thickness on the nucleation and seismicity is carefully analyzed in Section 4.

## 3. Terracotta Friction Law: Performance and Validation

The purpose of the current section is to demonstrate the capabilities of the proposed Terracotta friction law in capturing the complex frictional behaviour of clay-rich fault gouges. To this end, the analysis will make use of a variety of experiments on clay gouges conducted by Saffer and Marone (2003), Bedford et al. (2022), and Ashman and Faulkner (2023), which explore the shear response of gouge media under different normal stress values and controlled velocity steps, as well as the evolution of the gouge's thickness *h*.

As a start, we determine the corresponding material and system parameters of the Terracotta friction law. In particular, an initial guess is made based on previous values obtained by Wiebicke and Einav (2024), which are then iteratively updated to fit the experimental curves, as specified below for each series of tests. The initial conditions are set according to the available experimental values of the normal stress, $\sigma_n^0$, and gouge's thickness, $h^0$. In the absence of information regarding the thickness of potential shear band, the gouge thickness, $h^0$ is assumed to match the initial thickness of the tested samples. While selecting a smaller value for $h^0$ would result in higher shear strain rates and alter the rheological parameters, the conceptual conclusions should remain valid for the scenarios considered here. Furthermore, since in the relevant papers information was not available concerning the initial values of the solid fraction $\phi^0$, the following analysis will consider it as a fitting parameter against the experimental results. Finally, the elastic strain state is initialized by isotropically compressing the simulated sample from a reference stress-free configuration, see Equation 23. Table 1 presents the model parameters and initial conditions for the gouge materials considered hereinafter.

### 3.1. Experiments From Saffer and Marone (2003)

The tests conducted by Saffer and Marone (2003) consist of shearing of synthetic and natural gouges, in a double-direct shear geometry, at room humidity and temperature. Specimens were subjected to different levels of constant normal stress and sheared with varying velocities ranging from 1 to 200 µm/s. Focusing on clay-rich gouges, here we only study the tests on illite shale and smectite. The gouge layers were prepared by applying an initial compaction to the original sample resulting in an initial thickness $h^0$ equal to 2.25 and 2.75 mm for smectite and illite shale, respectively. Next, the gouge was sheared with frequent velocity stepping.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Shear response of illite shale and smectite gouges under constant normal stress $\sigma_n$ = 20 MPa, shown in terms of the friction coefficient, $\mu$ = $\tau/\sigma_n$, versus the shear strain, $\gamma$: (a) experiments (Saffer & Marone, 2003), and (b) model's predictions.</figcaption></figure>

The initial conditions for the Terracotta friction law consider the effective normal stress prescribed in the experiments, that is $\sigma_n$ = 20 MPa, the aforementioned values for $h^0$, while the solid fractions were chosen as $\phi^0$ = 0.57, 0.71 for the illite and smectite gouges, respectively.

Figure 4 displays the frictional behavior during the experiments and the Terracotta friction law given the parameters in Table 1 for the corresponding smectite and illite shale tests. Both the experiments and the predictions are characterized by an initially rapid increase in the friction coefficient during the loading phase, before peak friction coefficient at around $\gamma \approx 0.5 \div 1$. Post-peak the friction coefficient gradually decays towards a plateau, whose value depends on the gouge's material, exemplified by the critical state constant, *M*, and the shear velocity.

The Terracotta friction law shows a strong capacity in replicating the responses of both specimens, including not only the under velocity stepping during mature shearing, but also the peak and pre-peak values of the friction coefficient. The initial elastic response is considered to be entirely due to the elasticity of the gouge material since the stiffness of the loading apparatus is much larger (cf. Section 3.3 in Saffer & Marone, 2003). Despite the very rich and variegate velocity dependence of the friction coefficient observed in the experiments, the new friction law is able to predict and distinguish between the large and limited velocity strengthening behavior displayed by illite and smectite, respectively. The ability of the new law to adequately represent both pre- and post-peak behaviors moves beyond previously proposed rate-and-state and micromechanically inspired frictional laws, where the elasticity of the gouge is ignored.

Some discrepancies remain between the experiments and the simulations in the amplitude of the friction coefficient jumps. While the experiments display a predominantly symmetrical trend of the friction coefficient with respect to up- and down-steps of the slip velocity, *v*, the predictions show jumps of $\mu$ more pronounced along one direction (up-steps for illite, and down-steps for smectite). The reason lies in the evolution of the meso-temperature, with a quadratic decay back to its steady-state value, cf. Equation 18. This yields a corresponding quadratic dependence of the friction coefficient on the slip velocity, see Appendix B, making it difficult to accurately model the friction evolution under velocity steps of multiple orders of magnitude. Nevertheless, the predictions are found to be in good agreement with the experiments for the entire range of variation of the shear strain, in spite of the complex loading scenarios and the idealized boundary conditions (cf. direct shear).

To further illustrate the performance of the model under velocity-stepping conditions, Figure 5 shows a detailed comparison of the Terracotta predictions for experimental results on illite shale and smectite. Here, the gouges were
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Detailed view of velocity steps on illite shale (a–c) and smectite (d–f) under constant normal stress $\sigma_n$ = 100 MPa: (a, d) friction coefficient, $\mu$ (experiments after Saffer & Marone, 2003) and model's predictions, (b, e) gouge's thickness, $h$, and (c, f) meso-temperature, $T_m$ predicted by Terracotta friction law. The subscript and superscript $s$ in $\mu_s$, $h_s$, and $T_m^s$ denote values at steady-state.</figcaption></figure>

sheared under constant normal stress $\sigma_n$ = 100 MPa, and subjected to a tenfold increase of the velocity till a critical state has been achieved, followed by a tenfold decrease of the velocity (cf. Figures 6 and 7 in Saffer & Marone, 2003). Calibration of the material parameters for this particular protocol yields the following values: $c = 45 \cdot 10^6$ Ks, $a \equiv c$, $\eta = 53 \cdot 10^3$ K$^{-1}$s$^{-1}$, and $M = 0.86$ for illite, and $c = 11 \cdot 10^7$ Ks, $a \equiv c$, $\eta = 13 \cdot 10^3$ K$^{-1}$s$^{-1}$, and $M = 0.21$ for smectite, while all remaining parameters are kept as those listed in Table 1.

The Terracotta friction law successfully describes the velocity-strengthening behavior of the gouges, with only some minor differences in the predictions for smectite. Figure 5 also presents the evolution of the meso-temperature $T_m$ and the gouge thickness $h$. The meso-temperature responds almost instantaneously to changes in slip velocity, shifting toward a new stationary value and causing a peak in the friction coefficient before it relaxes to a new steady-state. A similar pattern is observed in the gouge thickness, which undergoes compaction under velocity down-steps and dilation under velocity up-steps.

## 3.2. Experiments From Bedford et al. (2022)

Bedford et al. (2022) performed experiments on fault gouge materials within a triaxial deformation apparatus in a direct-shear arrangement, involving both fully homogeneous and laterally heterogeneous gouges of kaolinite, quartz, and mixtures of them. The gouges, with an initial thickness $h^0 \approx 1$ mm, were subjected under a constant (effective) normal stress $\sigma_n$ = 40 MPa and sheared under relatively slow velocities between $0.3 \le \nu \le 3$ µm/s.

Here, we only focus on the clay-rich experiment that involved kaolinite, as presented in Figure 6. The initial thickness and normal stress are set to match the experimental values, while the initial solid fraction is chosen as $\phi^0 = 0.65$. This corresponds to an overconsolidation ratio of $p_c/p = 1$, where $p = \sigma_n$ at the initial state.

As for the previous tests, the evolution of the friction coefficient during the slip was provided, as well as the evolution of the gouge thickness. The latter was tracked by measuring the volume of fluid expelled from the sample by means of servo-control pumps, which revealed a substantial amount of consolidation. In particular, the inference of the thickness was made while assuming constant sliding area, such that the volumetric straining was
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison between experiments (Bedford et al., 2022) and the Terracotta friction law for kaolinite at normal stress $\sigma_n$ = 40 MPa. Evolution of (a) friction coefficient, $\mu$, (b) gouge thickness, $h$, and (c) normal velocity, $v_n$.</figcaption></figure>

approximated through changes only to the layer thickness without orthogonal straining. Importantly, the data of the thickness variations allows us to corroborate the volumetric nature of the Terracotta friction law.

The new friction law is calibrated based on the time evolution of the gouge thickness, $h$, presented in Figure 6b, before the application of velocity steps that initiates after slip displacement of 1.5 mm. The comparison of the model against the experimentally measured friction coefficient, gouge thickness, and normal velocity, $v_n$, is presented in Figure 6. Despite an underestimation of the peak value of the friction coefficient, the predictions beyond this peak are found to be in good agreement with the experimental results, and demonstrate capacity to capture the isotach phenomenon (Šuklje, 1969) at both up- and down-steps of the velocity. The evolution of the gouge thickness is well captured as well, thanks to the adopted calibration procedure. In particular, the model agrees with the experiments until the first downstep of the velocity, around 2 mm slip, where the friction law predicts a near constant sample thickness, in contrast with the slight continual compaction observed in the experiments. It is worth noticing that the intrinsic volumetric nature of the new friction law depends on the time evolution of two rates—slip velocity $v(t)$ as well as compaction velocity $v_n(t)$. To this purpose, we also inferred the experimental evolution of $v_n$, by means of a first-order finite difference scheme from the test data of the thickness, that is, $v_n(t) \approx v(t)\Delta h/\Delta d$, where $\Delta h$ and $\Delta d$ are the increments of the gouge's thickness and slip displacement. As shown in the corresponding Figure 6c, the predicted normal rate, $v_n$, well matches the experiments. Despite some minor differences mostly due to the resolution of the recorded evolution of the thickness, the Terracotta friction law correctly predicts the values of the compaction velocity before and during velocity stepping, with minor underestimation of the repeated upsteps of the shear velocity.

Finally, it is worth noticing that the minor differences between the predictions and the experiment may be due to the working assumption of material homogeneity in the current initial development. A further reason may be due to the idealised boundary conditions considered in the specification of the new friction law (see Section 2.3) as compared to the actual setup of the present experiment (Bedford et al., 2022), where soft silicone spacers were positioned at end of the gouge layers in order to accommodate longitudinal displacements without supporting any load.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** The variations in the friction coefficient, μ, of kaolinite in response to seven slip velocity steps under different constant normal stresses of σ<sub>n</sub> = 10, 25, 60 MPa, in terms of experiments (Ashman & Faulkner, 2023), and the Terracotta friction law. The friction coefficient curves are plotted at the same scale, around a mean value of μ = 0.33, with an offset of 0.15 along the y-axis for clarity.</figcaption></figure>

**3.3. Experiments From Ashman and Faulkner (2023)**

The experiments performed by Ashman and Faulkner (2023) use the same experimental setup of Bedford et al. (2022), again on kaolinite gouges in order to study their frictional behavior under velocity steps. However, in their study they additionally inspected the influence of the normal stress, σ<sub>n</sub>. Accordingly, we consider the exact same material and rheological parameters obtained from the calibration of the previous test, as listed in Table 1, and vary the initial solid fraction to keep a constant overconsolidation ratio of p<sub>c</sub>/p = 1.

Figure 7 compares the predictions provided by the Terracotta friction law against the experimental results. An overall good agreement is found, which is notable considering that no adjustment of the material parameters was made for the different normal stresses thanks to the dependence of the viscous component of the steady-state friction coefficient on σ<sub>n</sub> (cf. Appendix B). In particular, the decreasing amplitude of the isotachs at increasing normal stress is accurately predicted by the new friction law, with only minor differences visible under the least stressed configuration, that is, σ<sub>n</sub> = 10 MPa.

**3.4. Parametric Analysis**

This paragraph presents a parametric study that demonstrates the capacity of the Terracotta friction law in recovering a rich spectrum of rate-and-state responses. The setup mirrors the velocity-step protocol of Figure 5a for illite, under identical initial and boundary conditions. Here, the parametric study systematically explores the roles of (a) the rheological coefficients and (b) the elastic moduli, while focusing on their influence on the steady-state response and transient behavior during velocity jumps.

**3.4.1. Impact of the Rheological Coefficients**

To understand the phenomenological significance of the rheological constants, *c* (or *a*) and *η*, we perform two sets of tests, while holding all other parameters fixed. In the first set, Figures 8a–8c, we use a combination of constants *c* (= *a*) and *η* that maintains the same stationary viscous stress at a given strain rate by keeping *c* ∝ *η*¹/³. In the second set, Figures 8d–8f, we vary *c* and *η* while maintaining *c* ∝ *η*⁻¹ which affects the stationary viscous stress.

It is shown that by enforcing the proportionality c ∝ η¹/³, the stationary viscous component of the friction co-efficient is maintained constant during the initial loading phase, allowing any subsequent variations in the relaxation rate to be attributed exclusively to differences in the transient evolution of the meso-temperature. Figures 8a–8c demonstrate that higher values of *c* and *η* accelerate the relaxation of the transients of the meso-temperature and friction coefficient, while increasing the magnitude of compaction/dilation of the gouge layer, without perturbing the early time response.

Conversely, when *c* and *η* are varied while keeping the product *c*η fixed, the magnitude of the friction coefficient jumps increases with increasing *c* due to a greater component of the viscous stress, as shown in Figures 8d–8f.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The impact of the rheological parameters on the model response under velocity steps (cf. Figure 5). To this end, $c = a$ and $\eta$ are varied according to two protocols. The first considers $c^3/\eta = \bar{c}^3/\bar{\eta}$ with $\bar{c} = 45 \cdot 10^6$ Ks, $\bar{\eta} = 53 \cdot 10^3$ K$^{-1}$s$^{-1}$): (a) friction, (b) meso-temperature, and (c) gouge's thickness. The second maintains $c\eta = \bar{c}\bar{\eta}$ constant: (d) friction, (e) meso-temperature, and (f) gouge's thickness. All variables are shown as offsets with respect to the steady-state values preceding the first velocity jump for each simulation (refer to Figure 5).</figcaption></figure>

This in turn affects the underlying compaction and dilation of the gouge. The shape of the transient response to velocity jumps is not affected as the transient of the meso-temperature is kept identical during the loading.

These analyses demonstrate the ability of the Terracotta friction law to predict diverse rate-dependent phenomena (see also Wiebicke & Einav, 2024). The model predicts the characteristic behavior displayed by illite with rapid transitions of the stress between constant plateaus under step changes in the strain rate, a phenomenon often referred to as “isotach” behavior (Šuklje, 1969) in soil mechanics. The model also predicts the instantaneous jumps of the shear stress followed by a relaxation towards a new residual rate-dependent plateau as shown by smectite. Those rate-dependent phenomena are governed by one of the new state variables introduced in Terracotta, the meso-temperature. This is in contrast with classical rate-and-state frictional laws (Dieterich, 1979, 1981; Ruina, 1983) where the form of the rate-dependence is described in terms of (fixed) parameters (a and b). Here, the meso-temperature influences the range of predicted rate-dependent behaviors.

To better examine how $T_m$ controls the transient frictional response, Figure 9 presents magnified views of the evolution of $T_m$ and $\mu$, corresponding to the parametric study with $c \propto \eta^{1/3}$ that was shown in Figures 8a–8c. Although relative to the friction coefficient, $T_m$ rapidly attains its steady-state value (see Figure 8), the rate at which it evolves—governed by its evolution Equation 18—directly controls the time taken to establish the peak friction coefficient and consequently its value. For decreasing values of $c$, the transient dynamics of $T_m$ prolongs. As a result, we observe longer durations for the friction coefficient to establish its peak and increased peak values. Eventually, well after $T_m$ reaches its stationary value, the friction coefficient slowly relaxes due to the progressive reduction of the viscous stress.

### 3.4.2. Impact of the Elastic Coefficients

Next, we explore how the elasticity of the gouge through its elastic stiffness moduli, $\tilde{G}$ and $\tilde{K}$, influence the response. By maintaining the ratio $\tilde{G} \propto \tilde{K}$ fixed, we keep the critical stress ratios between $M_e$, $M_0$ and $M$ unchanged (Wiebicke & Einav, 2024), which ensures the requirements expressed in Equations 15 and 22 are met in
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The role of the meso-temperature in establishing the short-term, peak friction response. Here, the rheological parameters satisfy $c^3/\eta = \tilde{c}^3/\tilde{\eta}$, with different $c = a$ and $\eta$, as per Figure 8a: (left column) $c = 2\tilde{c}$, (center) $c = \tilde{c}$, and (right) $c = \tilde{c}/3$. Subfigures (a–c) refer to the velocity downstep and subfigures (d–f) display the velocity upstep.</figcaption></figure>

all the simulations. Figures 10a–10c displays the friction coefficient, the meso-temperature, and the gouge's thickness for both velocity-down and velocity-up steps. Panels (d–g) show zooms on the transient phases.

Increasing $\tilde{G}$ and $\tilde{K}$ results in larger magnitudes of the friction jumps and faster relaxations toward the new steady state. The physical elasticity of the material thus provides a new meaning to the extent of the so-called “critical slip distances” that are commonly lumped into a phenomenological constant in classical rate-and-state friction models (cf. $D_c$). By enforcing the proportionality $\tilde{G} \propto \tilde{K}$, the steady-state value and magnitude of the jumps of the meso-temperature remain unchanged. Instead, the gouge thickness evolution is strongly affected by the elasticity of the system: with higher jumps under increased elastic moduli $\tilde{G}$ and $\tilde{K}$. According to the present model, the extent of dilation (compaction) is a direct signature of the elasticity of the system, further suggesting a potential indirect measure of gouge elasticity in experiments.

## 4. Terracotta-Spring Model for Earthquakes Dynamics

This section extends the analysis to investigate the predictive capabilities of the Terracotta friction law in modelling earthquake nucleation, stick-slip instabilities, and their occurrence. We focus attention on exploring the influence of the physical parameters that are unique to this new friction law, such as slip thickness, meso-temperature, and solid fraction, and which cannot be effectively examined using conventional rate-and-state models (Dieterich, 1979, 1981; Kanamori & Brodsky, 2004; Pipping, 2019; Ruina, 1983). To this end, we connect the Terracotta friction law to a spring-dashpot representation of the host rock, as shown in Figure 1c. We then use this combined Terracotta-spring model to analyze the statistical dynamics of earthquakes and its energetical terms. In this model the mass, $m$, represents a mobilized mass sliding over a frictional gouge layer of thickness $h$ and length $L$, undergoing a slip displacement $d$ due to the far-field movement of the tectonic plates, which is prescribed by the slip velocity $v_\infty$. The mass is connected to a spring with stiffness $k = G_*L$, which represents the effective mobilized stiffness of the host rock, where $G_*$ is the corresponding shear modulus. This spring is connected in parallel to a dashpot with damping coefficient $\zeta$, equivalent to the apparent viscosity of the host rock. Subjected to an effective overburden normal stress $\sigma_n$, the frictional shear layer develops a shear stress $\tau$, governed by the evolution of Terracotta's state variables and the slip velocity, $v = \dot{d}$, see Section 2.3.

### 4.1. Governing Equations

Using force balance, the motion of the Terracotta-spring model is described as follows

$$m\dot{v} = k(d_\infty - d) + \zeta(v_\infty - v) - S\tau, \quad (27)$$
<figure><img src="image_12.png" /><figcaption>Figure 10. The impact of the elastic parameters on the model response under velocity steps (cf. Figure 5). To this end, $\tilde{G}$ and $\tilde{K}$ are varied by ensuring $\tilde{G}/\tilde{K} = \tilde{G}/\tilde{K}$, with $\tilde{G} = 400$ MPa and $\tilde{K} = 1040$ MPa: (a) friction, (b) meso-temperature, and (c) gouge's thickness. Subfigures (d) and (e) present magnified plots of the friction coefficient during the velocity lower and higher jumps, respectively. All variables are shown as offsets with respect to their respective steady-state values preceding the first velocity jump (a–e) and the second velocity jump (f, g).</figcaption></figure>

where $m = \rho_\ast L^3$ is the mobilised mass assuming a cuboid of length $L$, $S = L^2$ is the contact area, $d_\infty = v_\infty t$ is the displacement at the extremity of the spring-dashpot configuration, and $\zeta = 2\zeta\omega_n m$, where $\zeta$ is the damping ratio and $\omega_n = \sqrt{G_\ast/\rho_\ast}/L$ is the natural frequency of the system, assuming an equivalent linearly elastic isotropic medium for the host rock, while $\rho_\ast$ is its bulk density.

Next, we introduce the following dimensionless variables

$$ \hat{t} = \frac{t}{\mathcal{T}}, \quad \gamma = \frac{d}{D}, \quad \gamma_\infty = \frac{d_\infty}{D}, \quad \varepsilon_n = 1 - \frac{h}{D}, \qquad (28a) $$

$$ \gamma' = \frac{v}{v_\infty}, \quad \varepsilon_h = -\frac{v_n}{v_\infty}, \quad \Theta = \frac{T_m}{\Theta}, \qquad (28b) $$

where $D$, $\mathcal{T}$, and $\Theta$ are respectively the characteristic length, time, and meso-temperature, that is,

$$ D = h^0, \quad \mathcal{T} = \frac{D}{v_\infty}, \quad \Theta = \frac{1}{\mathcal{T}} \sqrt{\frac{c}{3\eta}}. \qquad (29) $$

Here, $D$ coincides with the initial value of the fault gouge thickness, $\mathcal{T}$ denotes the timescale over which shear strain accumulates within the gouge, and $\Theta$ corresponds to the ambient meso-temperature—that is, the steady-state value reached under a constant slip velocity $\mathcal{T}^{-1}$ and constant dilation/contraction of the fault gouge. Accordingly, the non-dimensional forms of Equations 27 and 24 read

$$ \gamma'' = \hat{k}(\gamma_\infty - \gamma) + \zeta(1 - \gamma') - \hat{N}\mu, \qquad (30a) $$
<table><caption>**Table 2**

Baseline Parameters for the Terracotta-Spring Model per Metre Unit Length of a Fault Gouge Segment</caption><thead><tr><th colspan="8">Material parameters</th><th colspan="4">In situ properties</th></tr><tr><th>K&#773; (GPa)</th><th>G&#773; (GPa)</th><th>M (-)</th><th>&omega; (-)</th><th>&lambda; (-)</th><th>&phi;<sub>l</sub> (-)</th><th>c (Ks)</th><th>&eta; (K<sup>-1</sup>s<sup>-1</sup>)</th><th>G<sub>*</sub> (GPa)</th><th>&rho;<sub>*</sub> (kg/m<sup>3</sup>)</th><th>&zeta; (-)</th><th>v<sub>&infin;</sub> (&mu;m/s)</th></tr></thead><tbody><tr><td>18.2</td><td>7</td><td>0.9</td><td>0.6</td><td>10</td><td>0.225</td><td>10<sup>5</sup></td><td>0.01</td><td>10</td><td>2.500</td><td>0 &divide; 1</td><td>10<sup>-5</sup> &divide; 50</td></tr><tr><th colspan="12">Initial conditions</th></tr><tr><th>&sigma;<sub>n</sub> (MPa)</th><td colspan="8"></td><th>&phi;<sup>0</sup> (-)</th><th colspan="2">h<sup>0</sup> (mm)</th></tr><tr><td>5 &divide; 640</td><td colspan="8"></td><td colspan="2">0.52 &divide; 0.9</td><td colspan="2">0.001 &divide; 10</td></tr></tbody></table>

Note. The rheological parameters are assumed to be equal, that is, $a \equiv c$. The far-field tectonic velocity, $v_{\infty}$, varies between 1 mm/year (~3 &middot; 10<sup>-5</sup> &mu;m/s) and 50 &mu;m/s.

$$ \mu = \frac{\hat{A}}{(1 - \varepsilon_n)^6} \varepsilon_v^e \gamma^e + \frac{\hat{B}}{1 - \varepsilon_n} \theta \gamma', \qquad (30b) $$

with

$$ \begin{aligned} \hat{k} &= \frac{\epsilon^2}{C}, & \hat{\epsilon} &= 2\zeta\sqrt{\hat{k}}, & \hat{N} &= \frac{\hat{k}}{\epsilon\kappa}, \\ \hat{A} &= \phi^{0.6} \frac{\hat{G}}{\sigma_n}, & \hat{B} &= \tau^{-1} \Theta \frac{2c}{3\Gamma \sigma_n}, & \hat{\epsilon} &= \frac{D}{L}, & \hat{C} &= \frac{\rho_* v_{\infty}^2}{G_*}, & \hat{\kappa} &= \frac{G_*}{\sigma_n}, \end{aligned} \qquad (31) $$

where $C$ is the Cauchy number representing the ratio of inertial and elastic forces through the host rock and $\epsilon$ is the thickness-to-length ratio of the fault gouge.

The dynamics of the Terracotta-spring model is governed by six non-dimensional variables—$\gamma$, $\epsilon_n$, $\gamma'$, $\epsilon'_n$, $\epsilon_v^e$, $\gamma^e$, and $\theta$—in terms of the non-dimensional time $\hat{l}$. The model's behavior is determined through the time integration of Equation 30, along with the evolution equations for the state variables $\epsilon_v^e$, $\gamma^e$, $\theta$, and $\delta_n$, as presented in Equation 25. We recall that the evolution of the normal displacement, $\delta_n$, is prescribed by ensuring a constant stress $\sigma_n$. The initial conditions for these state variables are given in Equation 23.

## 5. Terracotta-Spring Model: Results

The previous section has described the attachment of the Terracotta friction law with a spring-dashpot system representative of the host rock. Here we analyze the corresponding performance for modelling earthquake dynamics. Besides the material constants associated with the friction law in Section 2.3 and the initial conditions ($\sigma_n$, $h^0$, and $\phi^0$), five additional non-dimensional parameters of the spring-slider system must be specified: $T$, $\epsilon$, $C$, $\kappa$, and $\zeta$. Their values depend on the actual configuration of the fault gouge, either in laboratory experiments or in the field. To this end, in the following a parametric analysis is performed, as motivated by previous such explorations of rate-and-state friction laws (Y. Gu & Wong, 1991; Marone et al., 1990; Rice, 1985; Rice & Ruina, 1983; Rudnicki, 1988).

Here, we consider a clay-rich gouge with material parameters similar to those calibrated using direct shear tests (see Section 3) and by triaxial experiments by Wiebicke and Einav (2024). In particular, we select values in the range obtained from these calibrations (see Section 3 and Wiebicke & Einav, 2024), namely: 0.016 $\le$ $K$ $\le$ 26 GPa, 0.013 $\le$ $G$ $\le$ 13 GPa, 0.42 $\le$ $M$ $\le$ 1.3, 0.4 $\le$ $\omega$ $\le$ 0.65, 3.62 $\le$ $\lambda$ $\le$ 17.2, 0.05 $\le$ $\phi_l$ $\le$ 0.45, 7 $\le$ $c$ $\le$ 3 &middot; 10<sup>7</sup> Ks, and 0.01 $\le$ $\eta$ $\le$ 3 &middot; 10<sup>6</sup> K<sup>-1</sup>s<sup>-1</sup>. The material and in situ dimensional properties are given in Table 2, where the volumetric and deviatoric rheological constants are assumed equal, that is, $a \equiv c$. The following analysis is presented in terms of a unit fault gouge segment, $L$ = 1 m, and constant stiffness to focus attention on the new and key ingredients of the Terracotta-spring model. In most scenarios, the results are being presented in their non-
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Response of the Terracotta-spring model: (a) friction coefficient $\mu$ and non-dimensional slip velocity (or non-dimensional strain rate) $\gamma'$ as functions of the shear strain $\gamma$; (b) friction coefficient $\mu$ during the non-dimensional time $\hat{t}$ (b); (c) a focus view on the evolution of various components of the friction coefficient; and (d) a focus view highlighting the synchronicity of the meso-temperature $\theta$ and the non-dimensional thickness, $1 + \varepsilon_n$. Results refer to $v_\infty$ = 1 cm/year, and initial conditions: $h^0$ = 0.1 mm, $\phi^0$ = 0.71, and $\sigma_n$ = 40 MPa.</figcaption></figure>

dimensional form to further enable scaling of the relationships identified in this work. Stability analyses of the present model are expected to gather additional insights in the occurrence and magnitude of the predicted stick-slip instabilities (cf. Ciardo & Viesca, 2024; Y. Gu & Wong, 1991; J.-C. Gu et al., 1984; Rice & Ruina, 1983), but such investigation remains beyond the scope of the present work.

## 5.1. Benchmark

To demonstrate the predictive capabilities of the Terracotta-spring model, we consider a benchmark scenario where the parameters are set as follows: $v_\infty$ = 1 cm/year (equivalent to $\sim 3 \cdot 10^{-4}$ µm/s), $\zeta$ = 0 (undamped scenario), and the following initial conditions: $h^0$ = 0.1 mm, $\phi^0$ = 0.71, and $\sigma_n$ = 40 MPa. The initial solid fraction corresponds to an overconsolidation ratio equal to $p_c/p \sim 2.6$. Figure 11 presents the corresponding response of the Terracotta-spring analogue model, including the evolution of the friction coefficient, $\mu$, and non-dimensional slip velocity (or non-dimensional strain rate), $\gamma'$, as functions of the slip deformation, $\gamma$, and the non-dimensional time, $\hat{t}$, as well as the evolution of the meso-temperature, normal strain, slip, and friction coefficient.

Figures 11a and 11b showcase three key behaviors: (a) the initial elastic buildup of the friction coefficient, (b) the nucleation of earthquakes triggered by frictional instability with sustained stick-slip behavior, and (c) a transition to smooth and steady sliding. This sequence begins with an isolated seismic event followed by persistent and periodic stick-slip limit cycles. The first event is characterized by a stress drop of 7.4 MPa (with a friction drop $\Delta\mu \sim 0.18$) accompanied by a fault slip of 0.7 mm occurring in $\sim 30$ s, reaching a maximum slip velocity $v > 10^5 v_\infty$. The subsequent periodic events display stress drops of about 4 MPa, with 0.4 mm slip displacements developing over $\sim 20$ s. These periodic events have a return period of $\sim 16$ days. The stick-slip behaviour gradually converges to a steady, aseismic slip around $\hat{t} = 1000$, equivalent to 10 years, with the chosen non-dimensionalisation.
Examining Figures 11b–11d, the proposed friction law is effectively predicting the nucleation of earthquakes, during which the slip velocity increases until the inertia starts dominating and cause the Terracotta layer to accelerate to a peak velocity after which the friction drops dynamically. The friction drop can be decomposed into two contributions, elastic and viscous, which follows the relation $\tau = \tau^e + \tau^d$ based on the derivation in Section 2. The elastic component of the friction, $\mu^e = \tau^e/\sigma_n$, initially drops sharply, while the gouge compacts (decreasing $\varepsilon_n$), which leads to higher slip velocities and causes the viscous component, $\mu^d = \tau^d/\sigma_n$, to increase. The elastic component then rapidly rebuilds under constant normal strain, reaching a stationary value as the viscous component diminishes. At this juncture, the system appears to shift into an emergent stability, where the material compacts leading to an additional decrease of the friction. Throughout the stick-slip cycles, the model predicts persistent dilation-compaction cycles of the fault gouge, associated with increase-decrease of the shear zone thickness, $1 + \varepsilon_n$ (see Figure 11d). Additionally, the evolution of the meso-temperature, $\theta$, follows a similar pattern to the slip velocity, but also accounts for the fluctuating motion at the mesoscale related to the normal strain rate, $\varepsilon'_n$, and the energy flowing towards the microscale.

In summary, the proposed Terracotta-spring model enables to study the nucleation of earthquakes, with the development of well-defined and persistent stick-slip sliding, which is synchronised with the evolution of the meso-temperature and the evolution of the shear zone thickness, with an intermittent pattern. To investigate the transition from stick-slip cycle to smooth sliding, cf. Figure 11b, we considered a critically damped scenario ($\zeta = 1$). In this case, the model yields the exact same dynamic behavior, where persistent stick-slip cycles develop, with only a time-increasing lag in the recurrence of the periodic seismic events of the order of 0.03 at $\hat{t} = 100$. Therefore, in the following we shall consider the model without the dashpot using $\zeta = 0$, for simplicity.

## 5.2. Influence of Physical Parameters

Here, we investigate the influence of physical parameters and namely the slip zone thickness, the solid fraction, the normal stress, and the far-field velocity. While the influence of some of these parameters has been extensively studied using previous rate-and-state and micromechanically inspired friction laws, as well as using high-fidelity simulations including finite element and discrete element methods, the effects of other physical parameters remain challenging to capture with simplified one-dimensional models. The objective here is twofold: (a) to examine the impact of key physical quantities which were not well studied in the past, such as slip zone thickness, and (b) to assess whether the newly proposed friction law exhibits trends consistent with those observed in experiments but also in conventional rate-and-state friction laws, thereby validating its relevance in capturing known frictional mechanisms.

### 5.2.1. Slip Zone Thickness, *h*

In contrast with conventional state-of-the-art friction models that treat friction as a surface phenomenon, a key feature of the Terracotta friction law is the incorporation of a finite slip zone thickness. In the Terracotta-spring model, the thickness evolves during both the earthquake nucleation and slip, and its initial value appears to significantly influence the frictional instabilities, limit cycles, and stress drops associated with each seismic event. It is also worth noticing that the current model does not account for phenomena such as particle breakage, thermal pressurization, flash heating, or frictional melting, each of which is thought to influence the shear zone's thickness and fault surface roughness evolution (for an extensive review, see Collins-Craft et al., 2020; Faulkner et al., 2010; Stathas & Stefanou, 2023). Nonetheless, the results shed light on additional important phenomena that complement those listed above.

Following the experimental evidence collected and presented by Rice (2006), we select initial thicknesses between $10 \ \mu\text{m} \ge h^0 \ge 10 \ \text{mm}$. The parametric study is conducted per metre unit length of the fault gouge segment, leading to a non-dimensional thickness parameter $10^{-5} \ge \epsilon \ge 10^{-2}$. Figure 12 shows the phase portraits (left and center) of the friction coefficient with respect to the slip velocity, $\gamma'$, and the ratio of normal-to-slip velocity, $\varepsilon'_n/\gamma'$, as well as the friction coefficient as a function of the accumulated slip (right) for decreasing thicknesses (from top to bottom). For the relatively largest thickness, that is, $\epsilon = 10^{-2} (h^0 = 10 \ \text{mm})$, the model predicts an isolated seismic event followed by a smooth, steady sliding. When the initial thickness decreases to $\epsilon = 10^{-3} (h^0 = 1 \ \text{mm})$, the model exhibits stick-slip behaviour with rapidly decreasing amplitudes, which ultimately converge to the same fixed-point attractor as in the $\epsilon = 10^{-2}$ case. The attractor coincides with the steady-state
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Influence of the non-dimensional gouged shear zone thickness, $\epsilon$, on earthquake dynamics. Evolution of the friction coefficient with respect to the non-dimensional slip velocity (left), the normal-to-slip velocity (center), and shear strain (right) under decreasing shear zone thickness (from top to bottom): (a) $\epsilon = 10^{-2}$, (b) $\epsilon = 10^{-3}$, (c) $\epsilon = 10^{-4}$, and (d) $\epsilon = 10^{-5}$. Results refer to slip velocity of $v_\infty = 1$ cm/year and initial conditions of $\phi' = 0.71$ and $\sigma_n = 40$ MPa.</figcaption></figure>

friction (see Appendix B), which is approximately given by $\mu \sim M/\sqrt{3}$ owing to the relatively low slip velocities developed. For even smaller initial thicknesses, that is, $\epsilon \le 10^{-4}$ ($h^0 \le 100$ µm), the model predicts sustained stick-slip limit cycles with amplitudes inversely related to the thickness. In this latter scenario, the limit cycles also appear in the normal-to-slip velocity, demonstrating that the coupling between frictional slip and the underlying compaction/dilation of the fault gouge actively participates in the dynamics of earthquakes.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Influence of the non-dimensional gouged shear zone thickness, $\epsilon$, on earthquake characteristics. From left to right: (a) identification of the friction drops, $\Delta\mu_p$ and $\Delta\mu$, for stick-slip and smooth sliding, (b) evolution of the friction drops, and (c) return period, $\hat{t}_r$ ($t_r$), as a function of $\epsilon$. Results refer to slip velocity of $v_\infty$ = 1 cm/year and initial conditions of $\phi^0$ = 0.71 and $\sigma_n$ = 40 MPa.</figcaption></figure>

The results are further analyzed and summarized in Figure 13 which shows the evolution of the friction drops and their recurrence intervals as function of the gouged shear zone thickness. As this thickness decreases, the non-dimensional recurrence period of the seismic events $\hat{t}_r$ increases. However the physical (dimensional) recurrence interval, $t_r$, decreases at the decreasing of the thickness within the stick-slip region. The friction drops are separated into two categories: the primary, typically larger drop, $\Delta\mu_p$, and subsequent drops, $\Delta\mu$ that, in the case of stick-slip behavior, correspond to the limit cycle amplitude (also refer to Figure 11). The results suggest that there is a critical shear zone thickness that produces the largest stress drop. For larger $\epsilon$, the friction drops decrease rapidly as stick-slip transitions to steady, aseismic sliding. Similarly, as $\epsilon$ decreases, the friction drops for subsequent limit cycles approach zero, reflecting the decreasing amplitude of stick-slip events, while the primary drop reaches an asymptote, $\Delta\mu_p \sim 0.16$, independently of $\epsilon$.

It is worth noticing that predicted transition between steady slip and stick-and-slip behaviors is a direct result of the coupled evolution of the meso-temperature and the volumetric nature of the Terracotta friction law, under the imposed direct-shear boundary conditions. Different initial gouge's thicknesses affect the evolution of the meso-temperature and its steady-state value (cf. Appendix B) which, in turn, alters the evolution of the thickness and thus the solid fraction.

The above results gather new insights on how stick-slip behavior and nucleation of earthquakes are strongly related to the shear zone thickness. Note that the negative proportionality in the friction drops and return intervals at increasing initial thicknesses is in agreement with laboratory test of synthetic granular fault gouges (Lyu et al., 2019).

### 5.2.2. Solid Fraction, $\phi$

Next, we examine the impact of the initial solid fraction on the seismic response. Figure 14 presents these results in terms of the $\mu - \gamma'$ limit cycles, the magnitude of friction drops, and the return interval. Figure 14a shows that the limit cycle attractor remains independent of $\phi^0$, leading to a constant $\Delta\mu$ and return interval $\hat{t}_r$, respectively in Figures 14b and 14c. However, the initial solid fraction influences the primary friction drop, with two different trends emerging. For $\phi^0$ higher than the threshold for a normally consolidated gouge, $\phi_c$, the initial friction drop $\Delta\mu_p$ increases with increasing $\phi^0$, where $\phi_c \sim 0.65$ is obtained by requiring $p_c(\phi)/p = 1$. This aligns with the concept of critical state in soil mechanics, where denser materials display a more pronounced stress peak followed by strain softening before reaching critical state. In contrast, for looser states ($\phi^0 \le \phi_c$), the friction drop increases as $\phi^0$ decreases, due to stress adjustments during the elastic buildup that precedes the primary drop and enabling the compaction of the gouge.

### 5.2.3. Normal Stress, $\sigma_n$

The results in Figure 15 shows responses for various constant normal stress, $\sigma_n$ from 5 to 640 MPa. As $\sigma_n$ decreases, both the primary and cyclic friction drops increase, in agreement with experimental evidence (Numelin
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Influence of the solid fraction, $\phi^0$, on earthquake dynamics and characteristics under constant normal stress $\sigma_n = 40$ MPa: (a) phase portrait of the friction versus the non-dimensional slip velocity and limit cycle, (b) friction drops, and (c) return period in terms of $\phi^0$. Results refer to $v_\infty = 1$ cm/year and initial conditions of $h^0 = 0.1$ mm shear band thickness. The value of the solid fraction corresponding to the normally consolidated scenario is $\phi_c = 0.65$.</figcaption></figure>

et al., 2007; Saffer & Marone, 2003). In particular, Figure 15a shows that the limit cycles tend to expand toward higher cyclic friction values and contract at lower slip velocities as $\sigma_n$ decreases. The return period instead increases exponentially with increasing normal stress.

To further assess the predictive capabilities of the Terracotta-spring model relative to other models, simulations, and data from natural and induced seismicity, we temporarily move from the non-dimensional analysis to consider physical variations of dimensional parameters. Figure 16a shows an example of a slip-weakening curve obtained from the new model and definition of energetic parts of the budget per unit area during a seismic event, for $\sigma_n = 40$ MPa. As conventionally referred to in geophysics, the energy budget comprises the three contributions explained in Figure 16a and identified as: (a) released elastic energy, $E_R$, (b) fracture energy, $E_G$, and (c) frictional energy, $E_F$ (Kanamori & Brodsky, 2004; Kanamori & Rivera, 2006). For $\sigma_n = 40$ MPa, we find $E_G = 0.21$ J/m$^2$, $E_F = 170$ J/m$^2$, and $E_R = 0.04$ J/m$^2$.

Figure 16b shows a power law increase of the cyclic stress drops with normal stress, which well compares with previous discrete element method numerical simulations and experimental studies on gouges (Bai & Konietzky, 2022; Bedford & Faulkner, 2021; Passelègue et al., 2016). Here, we get $\Delta\tau \propto \sigma_n^{0.63}$. Furthermore, the model predictions for the relationships between the stress drops, $\Delta\tau$, and fracture energy, $E_G$, with respect to the critical slip, $\gamma_c$, align well with prior experimental studies conducted on smooth and gouge faults (Aubry et al., 2018; Ohnaka, 2003; Passelègue et al., 2016, 2017; Sellers et al., 2003), rock fracturing (Yoshimitsu et al., 2014), acoustic emission data (Goodfellow & Young, 2014), and discrete element method numerical simulations (Bai & Konietzky, 2022; Casas et al., 2023; Kato, 2012). Both $\Delta\tau$ and $E_G$ well match a power law with slip: $\Delta\tau \propto \gamma_c^{1.016}$ and $E_G \propto \gamma_c^{2.008}$. A first-order approximation of the fracture energy, given by the triangular area $(\tau_p - \tau_r)/\gamma_c$, suggests $E_G \propto \Delta\tau\gamma_c = \gamma_c^{2.016}$, which indicates minimal released energy, $E_R$.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Influence of the normal stress, $\sigma_n$, on earthquake dynamics and characteristics under a constant overconsolidation ratio $p/p_c \sim 4$: (a) stick-slip limit cycles, (b) friction drops, and (c) return period in terms of $\sigma_n$. Results refer to $v_\infty = 1$ cm/year, and initial conditions of $h^0 = 0.1$ mm.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Influence of the normal stress on fault energetics and statistics: (a) schematic representation of the energy budget; (b) stress drop as a function of the normal stress; and (c, d) fracture energy and stress drop over $\gamma_c$ defined in (a).</figcaption></figure>

### 5.2.4. Far-Field Velocity, $v_\infty$

Finally, we study the influence of the tectonic far-field velocity, $v_\infty$, by testing over six order of magnitude variations, that is, $10^{-5} \le v_\infty \le 50$ µm/s. Figure 17 summarizes the results which are in agreement with the well-known dependence of fault gouges on $v_\infty$ (Y. Gu & Wong, 1991, 1994; Kanamori & Brodsky, 2004; Scholz, 2019).

The far-field velocity determines the loading rate of the fault, and is found to directly influence the return period and indirectly affect the stress drops and stick-slip limit cycles. Higher velocities lead to faster stress accumulation, producing more frequent seismic events of lower magnitude. On the contrary, lower velocities yield slower loading rates, which in turns produce larger return periods, with seismic events characterised by larger stress drops. This is true both for the primary stress drop, $\Delta\mu_p$, and the persistent periodic drops, $\Delta\mu$. It is worth noticing that at vary large velocities, the stick-slip limit cycles tend to converge to a fixed-point attractor at $\mu \sim M\sqrt{3}$.

The return interval and the far-field velocity are found to follow a power law, that is, $v_\infty \propto t_r^{-1.08}$, as shown in Figure 18 and in agreement with existing natural seismicity data and laboratory tests (Beeler et al., 2001).

## 6. Conclusions

This study develops a new pathway to study the frictional behavior of fault gouges that is grounded on hydrodynamic principles from theoretical physics (Landau & Lifshitz, 2013) that considers the way energy flows in heterogeneous media through the use of thermal and meso-related temperatures (Einav & Liu, 2023; Jiang & Liu, 2009). In particular, we built upon the hydrodynamic Terracotta model of clay (Wiebicke & Einav, 2024). Initially formulated for triaxial loading, this paper established its tensorially generalized form, in order to solve the particular loading conditions of directly sheared fault gouges. Unlike previous models, this pathway ensures an objective evaluation that does not engineer the friction law against observational trends. As a result, the new friction law accounts for the volumetric response of the materials, and directly captures the dependence of fault gouge dynamics on thickness, normal stress and solid fraction, as well as explains for the first time analytically the buildups of an elastic stress towards first rupture.

The performance of the Terracotta friction law was validated against laboratory experiments of clay-rich fault gouges across various scenarios, including velocity stepping, varying normal stress, and gouge compaction (Ashman & Faulkner, 2023; Bedford et al., 2022; Saffer & Marone, 2003). By integrating this friction law with a conventional spring-dashpot representation of the host rock, the full seismic cycle was simulated, showcasing both shear stress accumulation and earthquake nucleation. The numerical results of this Terracotta-spring model successfully captured several well-established trends observed in laboratory experiments, advanced simulations, as well as natural and induced seismicity. Key result includes the model predictions of the finding by Lyu et al. (2019).

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Influence of the far-field velocity, $v_\infty$ on fault energetics and statistics: (a) stick-slip limit cycles, (b) friction drops, and (c) return period versus $v_\infty$. Results refer to initial conditions of $h^0 = 0.1$ mm, $\phi^0 = 0.71$, and $\sigma_n = 40$ MPa.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Influence of the far-field velocity, $v_\infty$: relationships between $v_\infty$ and return period, $t_r$.</figcaption></figure>

of the dependence of friction drops and return intervals on the initial gouge thickness. Furthermore, the same Terracotta-spring model offers insights into behaviors that await experimental studies, including: (a) the existence of a critical shear zone thickness that maximizes stress drops; and (b) the independence of the return periods on the initial solid fraction of the gouge.

The newly developed friction law represents a significant step beyond phenomenological rate-and-state laws, in that it enriches understanding of the physics of friction of fault gouges. Thanks to the physical generality of the hydrodynamic procedure (Jiang & Liu, 2009; Landau & Lifshitz, 2013) and the simplicity of Terracotta (Wiebicke & Einav, 2024), the use of the new friction law should be applicable to model general heterogeneous geomaterials, beyond just clay-rich gouges as done in this paper. As shown in Einav and Liu (2023) and Y. Chen et al. (2023), the hydrodynamic procedure can accurately describe multi-physics coupled phenomena, which have important implications on fault dynamics (see Barbot, 2022; Rattez et al., 2018; Rice, 2006; Stathas & Stefanou, 2023; Veveakis et al., 2012, among others).

While in this paper the initial thickness of the fault gouge was taken as a parameter, the evolution equation of the meso-temperature in the complete Terracotta model showcases dependence on its second gradient that naturally introduces an evolving length scale to continuum mechanics. With that, it is anticipated that the Terracotta constitutive model could help addressing both the genesis of faults—through the localization of strain, kinetic fluctuations (meso-temperature), and thermal pressurization (among others, see Stathas & Stefanou, 2023)—and the subsequent evolution of rate- and state-dependent earthquake cycles.

## Appendix A: Numerical Implementation

The numerical implementation of Terracotta friction law, Equations 24–26, is developed in Python, by solving the initial value problem specified by the evolution equations for the state variables and the time derivative of the shear stress, Equation 24, given the initial conditions in Equation 23, using the `solve_ivp` package from `SciPy` library (Virtanen et al., 2020). The time derivative of the shear stress is obtained by time differentiation of Equation 24:

$$ \dot{\tau} = \frac{\partial \tau}{\partial \varepsilon_v^e} \dot{\varepsilon}_v^e + \frac{\partial \tau}{\partial \dot{\gamma}^e} \dot{\gamma}^e + \frac{\partial \tau}{\partial T_m} \dot{T}_m + \frac{\partial \tau}{\partial h} v_n + \frac{\partial \tau}{\partial v} \dot{v}, \quad (A1) $$

where $\dot{\varepsilon}_v^e, \dot{\gamma}^e, \dot{T}_m$, and $v_n$ (or $\dot{h}$) are prescribed by Equations 25 and 26. The slip acceleration, $\dot{v}$, is approximated with a backward finite-difference scheme.

The solution of the initial value problem is driven by the prescribed time evolution of the shear strain rate, $\dot{\gamma}$ (or equivalently, the slip velocity, $v$) and the imposed value of the normal stress, $\sigma_n$. To guarantee a constant value of the normal stress, a root-finding algorithm is used to identify the normal-strain rate, $\dot{\varepsilon}_n$ (or, $v_n$), to yield

$$ \dot{\sigma}_n = \frac{\partial \sigma_n}{\partial \varepsilon_v^e} \dot{\varepsilon}_v^e + \frac{\partial \sigma_n}{\partial \varepsilon_s^e} \dot{\varepsilon}_s^e + \frac{\partial \sigma_n}{\partial \varepsilon_n^e} \dot{\varepsilon}_n^e + \frac{\partial \sigma_n}{\partial T_m} \dot{T}_m + \frac{\partial \sigma_n}{\partial h} v_n + \frac{\partial \sigma_n}{\partial v_n} \dot{v}_n = 0, \quad (A2) $$

where the rates of the state variables are functions of the state, as well as the shear and normal strain rates. Equation A2 is numerically solved using the `SciPy`'s `minimize` package.

Further notice that in all simulations, the numerical implementation of Equation 26b uses the initial thickness, $h^0$, rather than the current thickness, $h$, for the sake of simplicity. This assumption is confirmed by the results where the normal strain is found to be infinitesimal small, that is, $\lesssim 5\%$.

The system of ordinary differential equations defining the initial value problem is numerically integrated using a fifth-order implicit Runge-Kutta method from the Radau IIA family (Virtanen et al., 2020). Adaptive time
stepping is employed, with absolute error tolerances varying between 10⁻⁸ and 10⁻¹², to ensure that no significant differences are observed in the computed solutions.

The numerical implementation of the Terracotta-spring model, Equation 30, is analogous to that of the Terracotta friction law, with the only difference that the slip velocity, $\dot{\gamma}$, and acceleration $\ddot{\gamma}$, are defined by the solution of the balance equation of the Terracotta-spring-dashpot system, cf. Equation 30a. For more details, we refer to the code (Masi & Einav, 2024).

## Appendix B: Steady-State Friction Law

The steady-state limit of the Terracotta friction law can be found by imposing stationarity for the state variables under constant slip velocity, $v \ge 0$, and effective normal stress, $\sigma_n$. The stationarity of those variables is given when they Equation 25 vanish:

$$ \dot{\epsilon}_v^e = 0, \quad \dot{\gamma}^e = 0, \quad \dot{T}_m = 0, \quad \dot{h} = 0. \qquad (B1) $$

Under the above conditions the stationary values of the meso-temperature and gouge's thickness become

$$ T_m^s = \sqrt{\frac{c}{3\eta}} \frac{v}{h_s}, \quad h_s = \phi^0 h^0 \left( \frac{\omega p_l}{\sigma_n - T_m^s / \Gamma} \right)^{1/\lambda}, \qquad (B2) $$

where we considered most accurately and for all practical purposes $\sigma_n^e = p^e$ at steady-state, or equivalently that $\sigma_e^e = \sigma_n^e$.

Using the above relations, the steady-state friction stress reads

$$ \tau_s(T_m^s, h_s) = \frac{M}{\sqrt{3}} \left( \sigma_n - \frac{T_m^s}{\Gamma} \right) + B \sqrt{\frac{c}{3\eta}} \left( \frac{v}{h_s} \right)^2. \qquad (B3) $$

Furthermore, due to the second order dependence of $T_m^s$ on the strain rate in all considered scenarios, the thermodynamic pressure, $T_m^s/\Gamma$, is negligible compared to the elastic normal stress. Accordingly, we arrive at the following simple expression for the friction law at steady state:

$$ \tau_s(h_s) \approx \frac{M}{\sqrt{3}} \left( \sigma_n + \frac{B}{M} \sqrt{\frac{c}{\eta}} \left( \frac{v}{h_s} \right)^2 \right), \qquad (B4) $$

or equivalently for the steady-state friction coefficient

$$ \mu_s(h_s) \approx \frac{M}{\sqrt{3}} \left( 1 + \frac{B}{M\sigma_n} \sqrt{\frac{c}{\eta}} \left( \frac{v}{h_s} \right)^2 \right), \qquad (B5) $$

with $h_s \approx \phi^0 h^0 (\omega p_l / \sigma_n)^{1/\lambda}$.

The steady-state friction consists of an elastic contribution, proportional to the critical state line, $\propto M$, and a viscous term, proportional to the square of the velocity and inversely proportional to the applied normal stress, $\propto v^2/(h_s\sigma_n)$, where $h_s$ is the gouge thickness at steady-state and $\sigma_n$ is the normal stress. This quadratic dependence could be tuned in the future by altering the current linear dependence of the Onsager's rheological coefficients on $T_m$.
Notation

Einstein's summation is implied for repeated indices. i, j = 1, 2, 3

a, b, c volumetric, off-diagonal, and deviatoric plastic transport coefficients

d<sub>∞</sub>, d far-field and gouge's slip displacements

D characteristic length

e<sup>e</sup><sub>ij</sub> = ε<sup>e</sup><sub>ij</sub> - ε<sup>e</sup><sub>33</sub>δ<sub>ij</sub> deviatoric elastic strain tensor

f<sub>i</sub>, f<sup>m</sup><sub>i</sub> thermal and meso-related entropy currents

G̃, K̃ intrinsic shear and bulk moduli

G<sub>*</sub> host rock shear modulus

h gouged shear zone thickness

m mobilized mass

M critical stress ratio

L gouge's length

p = σ<sub>33</sub>/3 total pressure

p<sub>c</sub>(φ) isotropic compression line

p<sub>l</sub> constant defining the isotropic compression line

p<sup>d</sup> viscous pressure

p<sup>e</sup> elastic pressure

p<sup>T</sup> thermodynamic pressure

q = √3/2 s<sub>ij</sub>s<sub>ij</sub> deviatoric total stress (invariant)

q<sup>d</sup> viscous deviatoric stress

q<sup>e</sup> elastic deviatoric stress

R, R<sub>m</sub> thermal and meso-related entropy productions

s, s<sub>m</sub> thermal and meso-related entropies

s<sub>ij</sub> = σ<sub>ij</sub> - pδ<sub>ij</sub> deviatoric stress tensor

t time

T characteristic time

t̂ non-dimensional time

T, T<sub>m</sub> thermal temperature and meso-temperature

u, u<sup>e</sup>, u<sup>m</sup> total, elastic, and meso-related internal energy densities

v<sub>i</sub> velocity vector

v, v<sub>n</sub>, v<sub>∞</sub> slip, normal, and far-field velocities

X<sup>d</sup><sub>ijkl</sub>, X<sup>p</sup><sub>ijkl</sub> viscous and plastic transport coefficients

□<sup>0</sup> initial value of variable □

□̂ = ∂x/∂t derivative of variable □ with respect to time t

□' = ∂□/∂t derivative of the variable □ with respect to the non-dimensional time t̂
<table><tr><td>∇<sub>i</sub>□ = ∂□/∂x<sub>i</sub></td><td>gradient of the scalar variable □</td></tr><tr><td>∇<sub>i</sub>□<sub>i</sub> = ∂□<sub>i</sub>/∂x<sub>i</sub></td><td>divergence of the vectorial variable □<sub>i</sub></td></tr><tr><td>∇<sup>2</sup><sub>i</sub>□ = ∂<sup>2</sup>□/∂x<sub>i</sub><sup>2</sup></td><td>Laplacian of the scalar variable □</td></tr><tr><td>γ<sup>e</sup></td><td>elastic shear strain</td></tr><tr><td>γ̇, γ̇<sup>e</sup></td><td>total and elastic shear strain rates</td></tr><tr><td>Γ</td><td>fixed coefficient that sets the unit of the meso-temperature</td></tr><tr><td>η̄</td><td>energy sink coefficient</td></tr><tr><td>η</td><td>constant parameter specifying the energy sink</td></tr><tr><td>δ<sub>ij</sub></td><td>Kronecker delta</td></tr><tr><td>ε<sub>ij</sub><sup>e</sup></td><td>elastic strain tensor</td></tr><tr><td>ε<sub>v</sub><sup>e</sup> = ε<sub>ii</sub><sup>e</sup></td><td>volumetric stress invariant</td></tr><tr><td>ε<sub>s</sub><sup>e</sup> = √(2/3 ε<sub>ij</sub><sup>e</sup> ε<sub>ij</sub><sup>e</sup>)</td><td>deviatoric elastic strain invariant</td></tr><tr><td>ε̇<sub>ij</sub><sup>e</sup>, ε̇<sub>ij</sub><sup>p</sup>, ε̇<sub>ij</sub><sup>l</sup></td><td>elastic, total, and plastic strain-rate tensors</td></tr><tr><td>ε<sub>n</sub><sup>e</sup>, ε<sub>l</sub><sup>e</sup></td><td>normal and longitudinal elastic strains</td></tr><tr><td>ε<sub>n</sub></td><td>normal strain</td></tr><tr><td>ε̇<sub>n</sub>, ε̇<sub>l</sub></td><td>normal and longitudinal strain rates</td></tr><tr><td>μ</td><td>friction coefficient</td></tr><tr><td>φ</td><td>solid fraction</td></tr><tr><td>φ<sub>l</sub></td><td>solid fraction of the isotropic compression line, at p = 1 kPa</td></tr><tr><td>λ</td><td>slope of the isotropic compression line</td></tr><tr><td>ρ</td><td>bulk density</td></tr><tr><td>ρ<sub>s</sub></td><td>solid density</td></tr><tr><td>ρ<sub>*</sub></td><td>host rock's bulk density</td></tr><tr><td>σ<sub>ij</sub></td><td>total stress tensor</td></tr><tr><td>σ<sub>ij</sub><sup>d</sup></td><td>viscous stress tensor</td></tr><tr><td>σ<sub>ij</sub><sup>e</sup></td><td>elastic stress tensor</td></tr><tr><td>σ<sub>n</sub>, σ<sub>l</sub>, τ</td><td>normal, longitudinal, and shear stresses</td></tr><tr><td>Θ</td><td>characteristic meso-temperature</td></tr><tr><td>ω</td><td>position of the critical state line</td></tr><tr><td>ζ</td><td>damping ratio</td></tr></table>

## Data Availability Statement

The experimental data sets used for the validation of the model were obtained directly from the corresponding works (Ashman & Faulkner, 2023; Bedford et al., 2022; Saffer & Marone, 2003). The software for reproducing the results presented in this work is available at (Masi & Einav, 2024).
Acknowledgments

The authors would like to acknowledge the support of the Australia Research Council (ARC) under the Discovery Projects scheme (Grant agreement ID DP220101164: "Physics-informed hydrodynamic model for clay across scales"). Special thanks are also extended to Dr. Max Wiebicke for valuable discussions about the developed model. Open access publishing facilitated by The University of Sydney, as part of the Wiley - The University of Sydney agreement via the Council of Australian University Librarians.

References

Alaei, E., Marks, B., & Einav, I. (2021). A hydrodynamic-plastic formulation for modelling sand using a minimal set of parameters. *Journal of the Mechanics and Physics of Solids*, **151**, 104388. https://doi.org/10.1016/j.jmps.2021.104388

Ashman, I., & Faulkner, D. (2023). The effect of clay content on the dilatancy and frictional properties of fault gouge. *Journal of Geophysical Research: Solid Earth*, **128**(4), e2022JB025878. https://doi.org/10.1029/2022JB025878

Aubry, J., Passelègue, F. X., Deldicque, D., Girault, F., Marty, S., Lahfid, A., et al. (2018). Frictional heating processes and energy budget during laboratory earthquakes. *Geophysical Research Letters*, **45**(22), 12–274. https://doi.org/10.1029/2018GL079263

Bai, Q., & Konietzky, H. (2022). Numerical insight into energy partitioning during stick-slip events based on the framework of rate-and-state friction law. *Geophysical Research Letters*, **49**(9), e2021GL096677. https://doi.org/10.1029/2021GL096677

Barbot, S. (2019). Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and deterministic chaos in a single asperity fault. *Tectonophysics*, **768**, 228171. https://doi.org/10.1016/j.tecto.2019.228171

Barbot, S. (2022). A rate-state-and temperature-dependent friction law with competing healing mechanisms. *Journal of Geophysical Research: Solid Earth*, **127**(11), e2022JB025106. https://doi.org/10.1029/2022JB025106

Barras, F., & Brantut, N. (2025). Shear localisation controls the dynamics of earthquakes. *Nature Communications*, **16**(1), 711. https://doi.org/10.1038/s41467-024-55363-y

Bedford, J. D., & Faulkner, D. R. (2021). The role of grain size and effective normal stress on localization and the frictional stability of simulated quartz gouge. *Geophysical Research Letters*, **48**(7), e2020GL092023. https://doi.org/10.1029/2020GL092023

Bedford, J. D., Faulkner, D. R., & Lapusta, N. (2022). Fault rock heterogeneity can produce fault weakness and reduce fault stability. *Nature Communications*, **13**(1), 326. https://doi.org/10.1038/s41467-022-27998-2

Beeler, N., Hickman, S., & Wong, T.-F. (2001). Earthquake stress drop and laboratory-inferred interseismic strength recovery. *Journal of Geophysical Research*, **106**(B12), 30701–30713. https://doi.org/10.1029/2000JB900242

Bigaroni, N., Scuderi, M. M., Cappa, F., Guglielmi, Y., Nussbaum, C., Aldega, L., et al. (2023). Frictional properties of Opalinus Clay: Influence of humidity, normal stress and grain size on frictional stability. *Geophysical Journal International*, **233**(1), 211–228. https://doi.org/10.1093/gji/ggac457

Bouchbinder, E., & Langer, J. (2009). Nonequilibrium thermodynamics of driven amorphous materials. I. Internal degrees of freedom and volume deformation. *Physical Review E: Statistical, Nonlinear, and Soft Matter Physics*, **80**(3), 031131. https://doi.org/10.1103/PhysRevE.80.031131

Byerlee, J., & Summers, R. (1976). A note on the effect of fault gouge thickness on fault stability. *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*, **13**(1), 35–36. https://doi.org/10.1016/0148-9062(76)90226-6

Carpenter, B. M., Scuderi, M. M., Collettini, C., & Marone, C. (2014). Frictional heterogeneities on carbonate-bearing normal faults: Insights from the Monte Maggio fault, Italy. *Journal of Geophysical Research: Solid Earth*, **119**(12), 9062–9076. https://doi.org/10.1002/2014JB011337

Casas, N., Mollon, G., & Daouadji, A. (2023). Influence of grain-scale properties on localization patterns and slip weakening within dense granular fault gouges. *Journal of Geophysical Research: Solid Earth*, **128**(3), e2022JB025666. https://doi.org/10.1029/2022JB025666

Chen, J., & Niemeijer, A. R. (2017). Seismogenic potential of a gouge-filled fault and the criterion for its slip stability: Constraints from a microphysical model. *Journal of Geophysical Research: Solid Earth*, **122**(12), 9658–9688. https://doi.org/10.1002/2017JB014228

Chen, J., & Spiers, C. J. (2016). Rate and state frictional and healing behavior of carbonate fault gouge explained using microphysical model. *Journal of Geophysical Research: Solid Earth*, **121**(12), 8642–8665. https://doi.org/10.1002/2016JB013470

Chen, T., & Lapusta, N. (2009). Scaling of small repeating earthquakes explained by interaction of seismic and aseismic slip in a rate and state fault model. *Journal of Geophysical Research*, **114**(B1), B01311. https://doi.org/10.1029/2008JB005749

Chen, Y., & Einav, I. (2024). A hydrodynamic finite element model for chemo-mechanically loaded poroelastic materials. *Computers and Geotechnics*, **168**, 106161. https://doi.org/10.1016/j.compgeo.2024.106161

Chen, Y., Guillard, F., & Einav, I. (2023). A hydrodynamic model for chemo-mechanics of poroelastic materials. *Géotechnique*, **74**(13), 1–16. https://doi.org/10.1680/jgeot.22.00260

Ciardo, F., & Viesca, R. C. (2024). Non-linear stability analysis of slip in a single-degree-of-freedom elastic system with frictional evolution laws spanning aging to slip. arXiv preprint arXiv:2407.16846. https://doi.org/10.48550/arXiv.2407.16846

Collins-Craft, N. A., Stefanou, I., Sulem, J., & Einav, I. (2020). A cosserat breakage mechanics model for brittle granular media. *Journal of the Mechanics and Physics of Solids*, **141**, 103975. https://doi.org/10.1016/j.jmps.2020.103975

Daub, E. G., & Carlson, J. M. (2008). A constitutive model for fault gouge deformation in dynamic rupture simulations. *Journal of Geophysical Research*, **113**(B12), B12309. https://doi.org/10.1029/2007JB005377

Daub, E. G., Manning, M. L., & Carlson, J. M. (2010). Pulse-like, crack-like, and supershear earthquake ruptures with shear strain localization. *Journal of Geophysical Research*, **115**(B5), B05311. https://doi.org/10.1029/2009JB006388

De Gennes, P.-G., & Prost, J. (1993). *The physics of liquid crystals* (No. 83). Oxford University Press. https://doi.org/10.1063/1.2808028

Dieterich, J. H. (1979). Modeling of rock friction: I. Experimental results and constitutive equations. *Journal of Geophysical Research*, **84**(B5), 2161–2168. https://doi.org/10.1029/JB084iB05p02161

Dieterich, J. H. (1981). Constitutive properties of faults with simulated gouge. *American Geophysical Union*, **24**, 103–120. https://doi.org/10.1029/GM024p0103

Edwards, S. F., & Oakeshott, R. (1989). Theory of powders. *Physica A: Statistical Mechanics and its Applications*, **157**(3), 1080–1090. https://doi.org/10.1016/0378-4371(89)90034-4

Einav, I., & Liu, M. (2018). Hydrodynamic derivation of the work input to fully and partially saturated soils. *Journal of the Mechanics and Physics of Solids*, **110**, 205–217. https://doi.org/10.1016/j.jmps.2017.10.004

Einav, I., & Liu, M. (2023). Hydrodynamics of non-equilibrium soil water retention. *Water Resources Research*, **59**(1). https://doi.org/10.1029/2022WR033409

Evans, J. P. (1990). Thickness-displacement relationships for fault zones. *Journal of Structural Geology*, **12**(8), 1061–1065. https://doi.org/10.1016/0191-8141(90)90101-4

Falk, M. L., & Langer, J. S. (1998). Dynamics of viscoplastic deformation in amorphous solids. *Physical Review*, **57**(6), 7192–7205. https://doi.org/10.1103/PhysRevE.57.7192

Faulkner, D., Jackson, C., Lunn, R., Schlische, R., Shipton, Z., Wibberley, C., & Withjack, M. (2010). A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones. *Journal of Structural Geology*, **32**(11), 1557–1575. https://doi.org/10.1016/j.jsg.2010.06.009

Goldhirsch, I. (2008). Introduction to granular temperature. *Powder Technology*, **182**(2), 130–136. https://doi.org/10.1016/j.powtec.2007.12.002

Goodfellow, S., & Young, R. (2014). A laboratory acoustic emission experiment under in situ conditions. *Geophysical Research Letters*, **41**(10), 3422–3430. https://doi.org/10.1002/2014GL059965
Gu, J.-C., Rice, J. R., Ruina, A. L., & Simon, T. T. (1984). Slip motion and stability of a single degree of freedom elastic system with rate and state dependent friction. *Journal of the Mechanics and Physics of Solids*, **32**(3), 167–196. https://doi.org/10.1016/0022-5096(84)90007-3

Gu, Y., & Wong, T.-f. (1991). Effects of loading velocity, stiffness, and inertia on the dynamics of a single degree of freedom spring-slider system. *Journal of Geophysical Research*, **96**(B13), 21677–21691. https://doi.org/10.1029/91JB02271

Gu, Y., & Wong, T.-f. (1994). Nonlinear dynamics of the transition from stable sliding to cyclic stick-slip in rock. *Nonlinear dynamics and predictability of geophysical phenomena*, **83**, 15–35. https://doi.org/10.1029/GM083p0015

Haines, S. H., Kaproth, B., Marone, C., Saffer, D., & Van der Pluijm, B. (2013). Shear zones in clay-rich fault gouge: A laboratory study of fabric development and evolution. *Journal of Structural Geology*, **51**, 206–225. https://doi.org/10.1016/j.jsg.2013.01.002

Henann, D. L., & Kamrin, K. (2013). A predictive, size-dependent continuum model for dense granular flows. *Proceedings of the National Academy of Sciences*, **110**(17), 6730–6735. https://doi.org/10.1073/pnas.1219153110

Hetland, E. A., Simons, M., & Dunham, E. M. (2010). Post-seismic and interseismic fault creep I: Model description. *Geophysical Journal International*, **181**(1), 81–98. https://doi.org/10.1111/j.1365-246X.2010.04522.x

Ikari, M. J., Carpenter, B. M., & Marone, C. (2016). A microphysical interpretation of rate- and state-dependent friction for fault gouge. *Geochemistry, Geophysics, Geosystems*, **17**(5), 1660–1677. https://doi.org/10.1002/2016GC006286

Ikari, M. J., Niemeijer, A. R., & Marone, C. (2011). The role of fault zone fabric and lithification state on frictional strength, constitutive behavior, and deformation microstructure. *Journal of Geophysical Research*, **116**(B8), B08404. https://doi.org/10.1029/2011JB008264

Jiang, Y., & Liu, M. (2009). Granular solid hydrodynamics. *Granular Matter*, **11**(3), 139–156. https://doi.org/10.1007/s10035-009-0137-3

Jiang, Y., & Liu, M. (2017). Why granular media are thermal, and quite normal, after all. *The European Physical Journal E*, **40**, 1–16. https://doi.org/10.1140/epje/i2017-11497-4

Kamrin, K., & Koval, G. (2012). Nonlocal constitutive relation for steady granular flow. *Physical Review Letters*, **108**(17), 178301. https://doi.org/10.1103/PhysRevLett.108.178301

Kanamori, H., & Brodsky, E. E. (2004). The physics of earthquakes. *Reports on Progress in Physics*, **67**(8), 1429–1496. https://doi.org/10.1088/0034-4885/67/8/R03

Kanamori, H., & Rivera, L. (2006). Energy partitioning during an earthquake. *Geophysical Monograph Series*, **170**, 3–13. https://doi.org/10.1029/170GM03

Kato, N. (2012). Dependence of earthquake stress drop on critical slip-weakening distance. *Journal of Geophysical Research*, **117**(B1), B01301. https://doi.org/10.1029/2011JB008359

Landau, L. D., & Lifshitz, E. M. (2013). *Statistical physics* (Vol. 5). Elsevier. https://doi.org/10.1016/C2009-0-24487-4

Lapusta, N., Rice, J. R., Ben-Zion, Y., & Zheng, G. (2000). Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate- and state-dependent friction. *Journal of Geophysical Research*, **105**(B10), 23765–23789. https://doi.org/10.1029/2000JB900250

Lieou, C. K., & Langer, J. (2012). Nonequilibrium thermodynamics in sheared hard-sphere materials. *Physical Review E: Statistical, Nonlinear, and Soft Matter Physics*, **85**(6), 061308. https://doi.org/10.1103/PhysRevE.85.061308

Liu, M. (2021). Temperatures in grains and plasma. *Europhysics Letters*, **135**(3), 36003. https://doi.org/10.1209/0295-5075/ac1bc6

Lyu, Z., Rivière, J., Yang, Q., & Marone, C. (2019). On the mechanics of granular shear: The effect of normal stress and layer thickness on stick-slip properties. *Tectonophysics*, **763**, 86–99. https://doi.org/10.1016/j.tecto.2019.04.010

Makedonska, N., Sparks, D., Aharonov, E., & Goren, L. (2011). Friction versus dilation revisited: Insights from theoretical and numerical models. *Journal of Geophysical Research*, **116**(B9), B09302. https://doi.org/10.1029/2010JB008139

Maranic, Z., Guillard, F., Baker, J., Einav, I., & Marks, B. (2021). A granular thermometer. *Granular Matter*, **23**(2), 1–15. https://doi.org/10.1007/s10035-020-01081-3

Marone, C. (1998). The effect of loading rate on static friction and the rate of fault healing during the earthquake cycle. *Nature*, **391**(6662), 69–72. https://doi.org/10.1038/34157

Marone, C., Raleigh, C. B., & Scholz, C. H. (1990). Frictional behavior and constitutive modeling of simulated fault gouge. *Journal of Geophysical Research*, **95**(B5), 7007–7025. https://doi.org/10.1029/JB095iB05p07007

Marone, C., Vidale, J. E., & Ellsworth, W. L. (1995). Fault healing inferred from time dependent variations in source properties of repeating earthquakes. *Geophysical Research Letters*, **22**(22), 3095–3098. https://doi.org/10.1029/95GL03076

Masi, F., & Einav, I. (2024). Hydrodynamics of fault gouges from constitutive to physic of the friction: Terraquake [Software]. *Zenodo*. https://doi.org/10.5281/zenodo.14219386

Morgenstern, N., & Tchalenko, J. (1967). Microscopic structures in kaolin subjected to direct shear. *Géotechnique*, **17**(4), 309–328. https://doi.org/10.1680/geot.1967.17.4.309

Niemeijer, A. R., & Spiers, C. J. (2007). A microphysical model for strong velocity weakening in phyllosilicate-bearing fault gouges. *Journal of Geophysical Research*, **112**(B10), B10405. https://doi.org/10.1029/2007JB005008

Noda, H. (2008). Frictional constitutive law at intermediate slip rates accounting for flash heating and thermally activated slip process. *Journal of Geophysical Research*, **113**(B9), B09302. https://doi.org/10.1029/2007JB005406

Numelin, T., Marone, C., & Kirby, E. (2007). Frictional properties of natural fault gouge from a low-angle normal fault, Panamint Valley, California. *Tectonics*, **26**(2), TC2004. https://doi.org/10.1029/2005TC001916

Ohnaka, M. (2003). A constitutive scaling law and a unified comprehension for frictional slip failure, shear fracture of intact rock, and earthquake rupture. *Journal of Geophysical Research*, **108**(B2), 2080. https://doi.org/10.1029/2000JB000123

Onsager, L. (1931). Reciprocal relations in irreversible processes. i. *Physical Review*, **37**(4), 405–426. https://doi.org/10.1103/PhysRev.37.405

Passelègue, F. X., Latour, S., Schubnel, A., Nielsen, S., Bhat, H. S., & Madariaga, R. (2017). Influence of fault strength on precursory processes during laboratory earthquakes. In *Fault zone dynamic processes* (pp. 229–242). American Geophysical Union (AGU). https://doi.org/10.1002/9781119156895.ch12

Passelègue, F. X., Schubnel, A., Nielsen, S., Bhat, H. S., Deldicque, D., & Madariaga, R. (2016). Dynamic rupture processes inferred from laboratory microearthquakes. *Journal of Geophysical Research: Solid Earth*, **121**(6), 4343–4365. https://doi.org/10.1002/2015JB012694

Perfettini, H., & Molinari, A. (2017). A micromechanical model of rate and state friction: I. Static and dynamic sliding. *Journal of Geophysical Research: Solid Earth*, **122**(4), 2590–2637. https://doi.org/10.1002/2016JB013302

Pipping, E. (2019). Existence of long-time solutions to dynamic problems of viscoelasticity with rate-and-state friction. *ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik*, **99**(11), e201800263. https://doi.org/10.1002/zamm.201800263

Platt, J. D., Rudnicki, J. W., & Rice, J. R. (2014). Stability and localization of rapid shear in fluid-saturated fault gouge: 2. Localized zone width and strength evolution. *Journal of Geophysical Research: Solid Earth*, **119**(5), 4334–4359. https://doi.org/10.1002/2013JB010711
Pranger, C., Sanan, P., May, D. A., Le Pourhiet, L., & Gabriel, A.-A. (2022). Rate and state friction as a spatially regularized transient viscous flow law. *Journal of Geophysical Research: Solid Earth*, 127(6), e2021JB023511. https://doi.org/10.1029/2021JB023511

Rattez, H., Stefanou, I., Sulem, J., Veveakis, M., & Poulet, T. (2018). The importance of Thermo-Hydro-Mechanical couplings and microstructure to strain localization in 3D continua with application to seismic faults. Part II: Numerical implementation and post-bifurcation analysis. *Journal of the Mechanics and Physics of Solids*, 115, 1–29. https://doi.org/10.1016/j.jmps.2018.03.004

Rice, J. R. (1985). Constitutive relations for fault slip and earthquake instabilities. In *Instabilities in continuous media* (pp. 443–475). https://doi.org/10.1007/BF02590151

Rice, J. R. (2006). Heating and weakening of faults during earthquake slip. *Journal of Geophysical Research*, 111(B5), B05311. https://doi.org/10.1029/2005JB004006

Rice, J. R., & Ruina, A. L. (1983). Stability of steady frictional slipping. *Journal of Applied Mechanics*, 50(2), 343–349. https://doi.org/10.1115/1.3167042

Riley, D., Einav, I., & Guillard, F. (2023a). A consistent derivation of soil stiffness from elastic wave speeds. arXiv preprint arXiv:2312.01666. https://doi.org/10.48550/arXiv.2312.01666

Riley, D., Einav, I., & Guillard, F. (2023b). A constitutive model for porous media with recurring stress drops: From snow to foams and cereals. *International Journal of Solids and Structures*, 262, 112044. https://doi.org/10.1016/j.ijsolstr.2022.112044

Rudnicki, J. W. (1988). Physical models of earthquake instability and precursory processes. *Pure and Applied Geophysics*, 126(2–4), 531–554. https://doi.org/10.1007/BF00879009

Ruina, A. (1983). Slip instability and state variable friction laws. *Journal of Geophysical Research*, 88(B12), 10359–10370. https://doi.org/10.1029/JB088iB12p10359

Saffer, D. M., & Marone, C. (2003). Comparison of smectite-and illite-rich gouge frictional properties: Application to the updip limit of the seismogenic zone along subduction megathrusts. *Earth and Planetary Science Letters*, 215(1–2), 219–235. https://doi.org/10.1016/S0012-821X(03)00424-2

Scholz, C. H. (2019). *The mechanics of earthquakes and faulting*. Cambridge University Press. https://doi.org/10.1017/9781316681473

Scott, D. R., Marone, C. J., & Sammis, C. G. (1994). The apparent friction of granular gouge in sheared layers. *Journal of Geophysical Research*, 99(B4), 7231–7246. https://doi.org/10.1029/93JB03361

Scuderì, M. M., Carpenter, B. M., & Marone, C. (2014). Physicochemical processes of frictional healing: Effects of water on stick-slip stress drop and friction of granular fault gouge. *Journal of Geophysical Research: Solid Earth*, 119(5), 4090–4105. https://doi.org/10.1002/2013JB010641

Segall, P., & Rice, J. R. (1995). Dilatancy, compaction, and slip instability of a fluid-infiltrated fault. *Journal of Geophysical Research*, 100(B11), 22155–22171. https://doi.org/10.1029/95JB02403

Sellers, E. J., Kataka, M. O., & Linzer, L. M. (2003). Source parameters of acoustic emission events and scaling with mining-induced seismicity. *Journal of Geophysical Research*, 108(B9), 2418. https://doi.org/10.1029/2001JB000670

Stathas, A., & Stefanou, I. (2023). Fault friction under thermal pressurization during large seismic-slip: Numerical analyses and extension of the model of frictional slip. *International Journal of Mechanical Sciences*, 248, 108184. https://doi.org/10.1016/j.ijmecsci.2023.108184

Şuklje, L. (1969). *Rheological aspects of soil mechanics*. Wiley-Interscience.

Tse, S. T., & Rice, J. R. (1986). Crustal earthquake instability in relation to the depth variation of frictional slip properties. *Journal of Geophysical Research*, 91(B9), 9452–9472. https://doi.org/10.1029/JB091iB09p09452

Van den Ende, M., Chen, J., Ampuero, J.-P., & Niemeijer, A. (2018). A comparison between rate-and-state friction and microphysical models, based on numerical simulations of fault slip. *Tectonophysics*, 733, 273–295. https://doi.org/10.1016/j.tecto.2017.11.040

Veveakis, E., Sulem, J., & Stefanou, I. (2012). Modeling of fault gouges with cosserat continuum mechanics: Influence of thermal pressurization and chemical decomposition as coseismic weakening mechanisms. *Journal of Structural Geology*, 38, 254–264. https://doi.org/10.1016/j.jsg.2011.09.012

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., et al. (2020). SciPy 1.0: Fundamental algorithms for scientific computing in python. *Nature Methods*, 17(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2

Volpe, G., Collettini, C., Taddeucci, J., Marone, C., & Pozzi, G. (2024). Frictional instabilities in clay illuminate the origin of slow earthquakes. *Science Advances*, 10(26), eadn0869. https://doi.org/10.1126/sciadv.adn0869

Wiebicke, M., & Einav, I. (2024). A simple hydrodynamic model for clay. *Journal of the Mechanics and Physics of Solids*, 192, 105789. https://doi.org/10.1016/j.jmps.2024.105789

Wood, D. (1990). *Soil behaviour and critical state soil mechanics*. Cambridge University Press. https://doi.org/10.1017/CBO9781139878272

Yoshimitsu, N., Kawakata, H., & Takahashi, N. (2014). Magnitude-7 level earthquakes: A new lower limit of self-similarity in seismic scaling relationships. *Geophysical Research Letters*, 41(13), 4495–4502. https://doi.org/10.1002/2014GL060306

Yu, B., Chen, J., Spiers, C. J., Ma, S., Zhang, M., Qi, W., & Chen, H. (2024). Frictional properties of simulated fault gouges subject to normal stress oscillation and implications for induced seismicity. *Journal of Geophysical Research: Solid Earth*, 129(9), e2024JB029521. https://doi.org/10.1029/2024JB029521