Reproducing laboratory earthquakes with a discrete-continuum model

We present a novel numerical model allowing to tak e the best part of continuum-based and discrete modelling in a single framework. This mode l is applied to the reproduction of laboratory earthquakes in a high-pressure triaxial cell. It al lows to represent most of the relevant phenomena at stake, including elastic stress build-up during loading, f ast and slow sliding events, seismic waves emission in the surrounding elastic medium and evolution of fault g ou e on the sliding interface. We review here some illustrative results obtained with this model and p ropose further research avenues.


Introduction
Laboratory earthquakes have become a major tool for the understanding of seismic and aseismic sliding on natural faults [1, 2, and references therein]. Outstanding progress has been made in the instrumentation of such experiments, but for practical reasons, it remains difficult for experimentalists to acquire the full mechanical and kinematic processes at stake on faults, at high sampling rates in space and time. This is where numerical simulations can be helpful.
We present a numerical model which aims at reproducing triaxial experiments on saw-cut marble samples described in [1,2]. This model is in 2D and miniaturized when compared to its real-world counterpart, but includes relevant deformation mechanisms. The hanging wall and footwall of the sample are represented using a meshfree continuum approach, apart from the areas located in the neighbourhood (within a few hundreds of micrometers) of their contacting surfaces. In these areas, conforming polygonal grains bonded by a cohesive-zone model are implemented in a DEM framework [3]. Realistic boundary conditions (in terms of the elasticity of the loading system, of the absorption of the elastic waves and of the fluid pressure applied on the lateral boundaries) are introduced. Constitutive laws (both in the continuum and discrete parts) are calibrated based on experimental results found in the literature [4].

Simulations
As for the experiments reported in [1,2], the numerical model reproduces the deformation of a saw-cut sample (Fig. 1). The sample is 8 mm-high and 4 mm-wide (1/10 scale with respect to the reference experiment). Each part of the fault is continuous and elastic, excepting the first 200 µm of the contacting surfaces which are composed of discrete rigid and conforming polygons (typical size of 10 µm). They interact through a Cohesive Zone Model (CZM), which means that they reproduce well the elastic properties of the bulk rock but can also be damaged and separated from the surface if a threshold shear or tensile strength is overcome at their boundaries. The CZM parameters are calibrated using additional simulations reproducing published results of triaxial compression tests on marble [4] (Fig. 2).

Figure 2. Independent Calibration of the Cohesive Zone
Model parameters based on experiments reported in [4]. Low confining stress leads to localized and brittle failure while large confining stress promotes distributed damage The saw-cut sample is submitted to a lateral confining stress σ3 (set to 45, 90, or 180 MPa) and to vertical compression at constant speed, through two elastic blocks which are designed to mimic the deformability of the experimental apparatus and to dissipate elastic waves triggered by seismic events at the precut interface. Simulations are stopped for a vertical strain of 2%. The fault roughness is limited to the irregularities brought by the polygonal shapes of the rigid grains. The simulations are performed with the code MELODY [5], which allows mechanical coupling between large collections of rigid grains and continuous deformable bodies.

Stress history
As shown in Fig. 3, simulations are first characterized by the progressive increase of stress fields within the samples, combined with damage on the first grain layers at the interface. At a critical differential stress, sliding events occur and stress values drop and become more heterogeneous. Time-series of stress and strain tensors show typical laboratory earthquake cycles, illuminating successive accumulation and release of shear stress on the fault (Fig. 4). The case σ3=45 MPa leads to a large number of low-amplitude stress drops, while the case σ3=180 MPa is characterized by a smaller number of large stress drops. The case σ3=90 MPa shows an intermediate behaviour. When expressed in terms of friction coefficient on the fault (Fig. 5)

Sliding velocity history
Fault sliding is recorded for each simulation at a high sampling rate, both in space and time. This is performed by monitoring displacements and velocities of the surface nodes of the continuum parts of the two half blocks and by projecting them on the mean fault plane. It provides a detailed velocity map of the fault for the three confining pressures (Fig. 6). A large number of events can be observed and characterized. Dedicated procedures are used for detection and characterization of each sliding event. To illustrate this database, we present here three characteristic events of increasing size. All of them occurred under a confining stress σ3=180 MPa. The first illustrative event (Fig. 7a) is of small amplitude, with a stress drop of ~0.4 MPa and a sliding distance of ~0.3 µm. This event originates in the middle part of the fault and propagates towards its extremities.  This results in the progressive creation of a third-body (i.e. granular gouge) layer separating the two walls of the fault.

Fault interface damage evolution
This view is confirmed by damage profiles across the fault, as provided in Fig. 8 (lower part). We observe that the central part of the fault is constituted of completely damaged material. Past this gouge layer, damage progressively decreases from complete degradation (damage=1) to almost intact rock (damage close to 0.2) following a decreasing exponential pattern. The influence of the confining pressure is highlighted. For values of σ3 increasing from 45 to 180 MPa, the thickness of the gouge layer evolves from ≈40 µm to ≈70 µm, respectively. The decrease of damage with the distance from the fault is also faster for low confining stresses than for larger ones.

Perspectives
The model described in this study is able to reproduce typical seismic cycles (frequency and amplitude of events) only relying on independent calibrations of its parameters, and so without using any ad-hoc slipweakening or rate-and-state friction law at contacts between elements composing the fault. Upon further use and interpretations, it is thus expected to provide valuable information on the local phenomena at stake during laboratory earthquakes.
Future work will be dedicated to overcome the current limitations of the model by increasing its size, by introducing a realistic roughness at the fault surface, and by implementing plasticity in the bulk rock continuum model, observed in [1,2]. More data about stress redistributions and fault damage evolution, as well as detailed statistics on acoustic and slip events, will also be gathered in order to analyse more deeply the influence of gouge emission on the fault rheology.