Extracting microphysical fault friction parameters from laboratory and field injection experiments

Human subsurface activities induce significant hazard by (re-)activating slip on faults, which are ubiquitous in geological reservoirs. Laboratory and field (decametric-scale) fluid injection experiments provide insights into the response of faults subjected to fluid pressure perturbations, but assessing the long-term stability of fault slip remains challenging. Numerical models offer means to investigate a range of fluid injection scenarios and fault zone complexities and require frictional parameters (and their uncertainties) constrained by experiments as an input. In this contribution, we propose a robust approach to extract relevant microphysical parameters that govern the deformation behaviour of laboratory samples. We apply this Bayesian approach to the fluid injection experiment of Cappa et al. (2019) and examine the uncertainties and trade-offs between parameters. We then continue to analyse the field injection experiment reported by Cappa et al. (2019), from which we conclude that the fault-normal displacement is much larger than expected from the adopted microphysical model (the Chen–Niemeijer–Spiers model), indicating that fault structure and poro-elastic effects dominate the observed signal. This demonstrates the importance of using a microphysical model with physically meaningful constitutive parameters, as it clearly delineates scenarios where additional mechanisms need to be considered.

1 Introduction

Induced seismicity is of primary concern in human subsurface activities, including geothermal energy production, wastewater and CO₂ injection, and hydrocarbon extraction (Ellsworth, 2013). Seismicity triggered around injection sites is generally attributed to elevated pore fluid pressures, which lower the clamping stress that keeps the fault locked (Elsworth et al., 2016). Additionally, recent field injection tests at a decametric scale reveal the importance of aseismic creep in driving seismicity (Duboeuf et al., 2017), and long-range poro-elastic effects and earthquake interactions have been inferred to trigger seismicity well beyond the extent of the stimulated region (Catalli et al., 2016; Goebel and Brodsky, 2018; Schoenball and Ellsworth, 2017). To better assess the earthquake hazard associated with the injection and extraction of geo-fluids, potential mechanisms underlying the nucleation of induced seismic events need to be identified.

Laboratory experiments provide the means to investigate the mechanisms for (unstable) fault slip at high resolution under well-controlled conditions (e.g. Kaproth et al., 2016; Scuderi et al., 2016, 2017; Tenthorey et al., 2003). Many laboratory studies report their results in terms of rate-and-state friction (RSF; Dieterich, 1979; Ruina, 1983) parameters, which may serve as input for numerical modelling studies (Cubas et al., 2015; Kroll et al., 2017; McClure and Horne, 2011; Noda et al., 2017). Unfortunately, it is well-established that the RSF parameters depend on a plethora of thermodynamic conditions (Blanpied et al., 1998; Boulton et al., 2019; Chester, 1994; He et al., 2016; Hunfeld et al., 2017), including fluid pressure (Cappa et al., 2019; Sawai et al., 2016; Scuderi et al., 2016), which needs to be accounted for when attempting to extrapolate laboratory measurements to nature through RSF-based numerical models. The relationships between RSF parameters and observable quantities (such as porosity, grain size, or fluid chemistry) are not well understood, and so great care must be taken when generalising laboratory results to natural systems.

As an alternative approach, decametric-scale fluid injection tests allow one to probe the response of a tectonic fault to fluid pressure perturbations under in situ conditions (Derode et al., 2015; Duboeuf et al., 2017; Guglielmi et al., 2015; Rivet et al., 2016). While these tests provide more direct insights into the (potentially seismic) behaviour of the fault, they are also more complicated to interpret owing to the complexity inherent to natural faults. Generalisation of the results and extrapolation to other fault or reservoir conditions is therefore challenging. Moreover, fluid injection rates and volumes are limited by regulatory restrictions, which inhibits a comparison with systems characterised by larger injection volumes and rates. Numerical models remain essential to investigate faults in this context (e.g. Dempsey and Riffault, 2019; Rutqvist et al., 2007; Wynants-Morel et al., 2020), which in turn rely on constraints offered by laboratory experiments.

In the present study, we reinterpret the laboratory and decametric-scale fluid injection experiments reported by Cappa et al. (2019) in the framework of the Chen–Niemeijer–Spiers (CNS) microphysical model (Chen and Spiers, 2016; Niemeijer and Spiers, 2007). To this end, we propose a robust approach for the extraction of the CNS microphysical parameters from laboratory or field observations based on the relation between fault dilatancy and shear slip and the temporal evolution of the slip rate. In this Bayesian approach, we examine the uncertainties associated with each parameter, and the trade-offs between parameters, which are both important for choosing suitable parameter ranges for numerical modelling efforts. Lastly, we discuss the limitations of and perspectives offered by the adopted microphysical model in the context of induced seismicity modelling.

2 Methods

2.1 The Chen–Niemeijer–Spiers model

To describe the observed laboratory observations of Cappa et al. (2019) in terms of micro-physical quantities, we adopt the Chen–Niemeijer–Spiers (CNS) model proposed by Niemeijer and Spiers (2007) and extended by Chen and Spiers (2016). In the following section, we briefly summarise the basic mechanics of this microphysical model and the numerical implementation adopted in this study. For a detailed derivation and discussion of this model, we refer to the original works of Niemeijer and Spiers (2007) and Chen and Spiers (2016) (see also Verberne et al., 2020, this special issue).

Firstly, the CNS model considers a representative elementary volume of fault gouge of thickness L and porosity ϕ, which is subjected to an effective normal stress σ_e (i.e. total normal stress minus the fluid pressure) and shear stress τ. In response to this state of stress, the gouge deforms internally through parallel operation of dilatant granular flow and one or more non-dilatant creep mechanisms. The timescales considered in the present study are too short (of the order of seconds to minutes) to justify a detailed consideration of the non-dilatant creep component, and hence we focus purely on the granular flow component. As will be shown later, this simplification is well-warranted by the laboratory observations. In line with this assumption, the shear and volumetric deformation of the fault gouge can be described as

$$\begin{array}{}\text{(1a)}& {\displaystyle \frac{\mathrm{d}\mathit{\delta }}{\mathrm{d}t}}& {\displaystyle }=V=L{\dot{\mathit{\gamma }}}_{\text{gr}},\text{(1b)}& {\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}t}}& {\displaystyle }=-\left(\mathrm{1}-\mathit{\varphi }\right){\dot{\mathit{\epsilon }}}_{\text{gr}}=\tan \mathit{\psi }\left(\mathrm{1}-\mathit{\varphi }\right){\dot{\mathit{\gamma }}}_{\text{gr}}.\end{array}$$

Here, V denotes the rate of slip on the fault δ and ${\dot{\mathit{\gamma }}}_{\text{gr}}$ and ${\dot{\mathit{\epsilon }}}_{\text{gr}}$ are the shear and volumetric strain rate of granular flow, respectively (compression defined positive). We consider only fault-normal volumetric strains (i.e. no fault-parallel expansion/contraction). The amount of volumetric deformation associated with an increment of shear strain is described by the dilatancy angle tan ψ, i.e. $\mathrm{d}{\mathit{\epsilon }}_{\text{gr}}=-\tan \mathit{\psi }\mathrm{d}{\mathit{\gamma }}_{\text{gr}}$, and is given by (Niemeijer and Spiers, 2007)

$$\begin{array}{}\text{(2)}& \tan \mathit{\psi }=\mathrm{2}H\left({\mathit{\varphi }}_{\mathrm{c}}-\mathit{\varphi }\right),\end{array}$$

where H is a geometric constant of order 1 and ϕ_c is referred to as the “critical-state” porosity, i.e. the maximum attainable porosity of the gouge. The parameter H represents how much dilatancy is involved when grains are sliding past one another and is likely affected by grain shape, angularity, and size distribution. Based on a first-order geometric analysis, Niemeijer and Spiers (2007) estimated that the maximum dilatancy angle at zero porosity is $\tan \mathit{\psi }=\sqrt{\mathrm{3}}$, which puts an upper bound on $H<\sqrt{\mathrm{3}}/\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}$. Likewise, the critical-state porosity ϕ_c is likely not a universal constant. Nonetheless, in the absence of tight theoretical constraints on H and ϕ_c, we treat these quantities as constant parameters.

The rate of granular flow is itself a function of stress and porosity and can be written as (Chen and Spiers, 2016)

$$\begin{array}{}\text{(3)}& {\dot{\mathit{\gamma }}}_{\text{gr}}={\dot{\mathit{\gamma }}}_{\text{gr}}^{*}\exp \left({\displaystyle \frac{\mathit{\tau }\left[\mathrm{1}-{\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{*}\tan \mathit{\psi }\right]-{\mathit{\sigma }}_{\mathrm{e}}\left[{\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{*}+\tan \mathit{\psi }\right]}{\stackrel{\mathrm{̃}}{a}\left[{\mathit{\sigma }}_{\mathrm{e}}+\mathit{\tau }\tan \mathit{\psi }\right]}}\right).\end{array}$$

The reference grain boundary friction coefficient ${\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{*}$ corresponds with a shear strain rate ${\dot{\mathit{\gamma }}}_{\text{gr}}^{*}$, and $\stackrel{\mathrm{̃}}{a}$ is a proportionality constant for the logarithmic velocity dependence of the grain boundary friction $\stackrel{\mathrm{̃}}{\mathit{\mu }}$, given by

$$\begin{array}{}\text{(4)}& \stackrel{\mathrm{̃}}{\mathit{\mu }}={\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{*}+\stackrel{\mathrm{̃}}{a}\ln \left({\displaystyle \frac{{\dot{\mathit{\gamma }}}_{\text{gr}}}{{\dot{\mathit{\gamma }}}_{\text{gr}}^{*}}}\right).\end{array}$$

We highlight that ${\dot{\mathit{\gamma }}}_{\text{gr}}$ is exponentially sensitive to the fluid pressure p through the effective stress ${\mathit{\sigma }}_{\mathrm{e}}=\mathit{\sigma }-p$, and so the CNS model predicts an acceleration of V upon an increase in the fluid pressure. Moreover, the experiments analysed in this study are conducted at constant shear stress, so that a force balance (which typically takes the place of Eq. 1a) is not required.

In the present study, we treat the laboratory sample as a single degree-of-freedom (spring block) system, with uniform porosity and an internal state of stress. This implies that the fluid pressure is considered to be uniform and constant throughout the sample, with no coupling between volumetric deformation and fluid pressure. This assumption is valid for samples with sufficiently high permeability, such that the characteristic timescale of fluid diffusion is smaller than the timescale of deformation. In other words, the sample is assumed to be in equilibrium with the externally applied fluid pressure (“drained”) at all times. In the laboratory experiments of Cappa et al. (2019), the gouge permeability was estimated to be above the intrinsic permeability of the apparatus (10⁻¹⁴ m²), so the sample can be considered to be drained. For low-permeability gouges, such as shales (Scuderi and Collettini, 2018), coupling between volumetric deformation and fluid pressure needs to be considered (e.g. Segall and Rice, 1995).

2.2 Microphysical parameter inversion procedure

In the simplified CNS framework laid out above, the dynamics of the system are fully governed by L, H, ϕ_c, $\stackrel{\mathrm{̃}}{a}$, and ${\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{*}$ (which simultaneously constrains ${\dot{\mathit{\gamma }}}_{\text{gr}}^{*}$), for a given state of stress and initial porosity. In principle, the forward model given by Eqs. (1a) and (1b) can be solved iteratively and used to invert laboratory measurements for these constitutive parameters. However, owing to the exponential sensitivity of V to ϕ through ${\dot{\mathit{\gamma }}}_{\text{gr}}$, such inversion procedure is unstable and ill-posed. As an alternative, we propose a two-step inversion procedure that robustly constrains the constitutive parameters. Firstly, we rewrite Eq. (1b) as

$$\begin{array}{}\text{(5)}& \mathrm{d}\mathit{\varphi }={\displaystyle \frac{\mathrm{2}H}{L}}\left({\mathit{\varphi }}_{\mathrm{c}}-\mathit{\varphi }\right)\left(\mathrm{1}-\mathit{\varphi }\right)\mathrm{d}\mathit{\delta },\end{array}$$

where dδ=Vdt is an increment of slip across the fault. While we recognise that L varies with ϕ, integration of Eq. (5) does not yield an analytical solution when taking L=f(ϕ). Fortunately, as will be shown later, we find that the inferred variations in L are of the order of 10 %–20 % of the absolute value of L, warranting a first-order approximation of a constant value of L. By integrating the above relation from the initial porosity ϕ₀ up to ϕ (cf. van den Ende et al., 2018) and recognising that $\mathrm{\Delta }L/L=\left(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}}\right)/\left(\mathrm{1}-\mathit{\varphi }\right)$ (for L≈L₀), we obtain an expression for the dilatancy ΔL as a function of slip δ:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{\Delta }L}{L}}={\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{c}}-{\mathit{\varphi }}_{\mathrm{0}}}{\mathrm{1}-{\mathit{\varphi }}_{\mathrm{c}}}}\left[\mathrm{1}-\exp \left(-\mathrm{2}H{\displaystyle \frac{\mathit{\delta }}{L}}\left[\mathrm{1}-{\mathit{\varphi }}_{\mathrm{c}}\right]\right)\right].\end{array}$$

This expression already provides sufficient means to constrain the constitutive parameters L, H, ϕ_c, and the initial condition ϕ₀ without numerically solving the full forward model given by Eqs. (1a) and (1b). The second step of the inversion involves constraining the remaining parameters $\stackrel{\mathrm{̃}}{a}$ and ${\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{*}$ by comparing Eq. (1a) with the laboratory measured slip rate. Since the slip rate can span orders of magnitude, we perform the inversion in terms of $\ln \left({\dot{\mathit{\gamma }}}_{\text{gr}}\right)$ (and correspondingly ln (V) measured during the experiment), which renders a more stable inversion task.

Since the proposed inversion protocol does not involve numerically solving a forward model, a single evaluation of either Eqs. (6) or (3) yields a sample of the posterior distribution, hence permitting extensive random sampling. To inspect the trade-offs between parameter values and their uncertainties, we cast the protocol above in a Bayesian inversion procedure, in which we estimate the posterior distributions $P\left(m=\left\{L,H,{\mathit{\varphi }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{c}}\right\}|\mathit{\delta },\mathrm{\Delta }L\right)$ and $P\left(\stackrel{\mathrm{̃}}{a},{\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{*}|\mathit{\delta },V,{\mathit{\sigma }}_{\mathrm{e}},\mathit{\tau },m\right)$ separately. We assume a uniform prior distribution over a bounded range of admissible parameter values, and a Gaussian likelihood with an unknown data variance ν² that is simply treated as a nuisance parameter and co-inverted. The posterior distributions are sampled using an affine invariant Markov chain Monte Carlo ensemble sampler as implemented in the Python emcee package (Foreman-Mackey et al., 2013). While it is also possible to estimate the posterior distributions from numerically solving the forward problem, each forward model evaluation from t=0 up to the point where V>1 mm s⁻¹ takes several tens of seconds on a single CPU. The practical reason for this is that the fault is critically stressed and hence requires small time step evaluations to ensure sufficient numerical accuracy and stability.

3 Analysis of fluid injection tests of Cappa et al. (2019)

3.1 Laboratory experiment

We apply the above procedure to the laboratory fluid injection experiment performed by Cappa et al. (2019) – see Fig. 1. In this experiment, a carbonate gouge sample was subjected to a constant shear stress of τ=1.2 MPa and a total normal stress of σ=5 MPa. The fluid pressure was increased step-wise every 150 s with steps of 0.5 MPa, until the sample “failed” macroscopically at a fluid pressure of p=3.5 MPa. Prior to the final stage of pressurisation, only negligible amounts of slip were measured, and hence we focus our inversion efforts on the final stage of the experiment in which the sample measurably accelerated. Additionally, through the stage of fluid injection, no gouge compaction was measured, supporting our assumption made prior to Eqs. (1a) and (1b) that the time-dependent creep rate is negligible compared to the rate of granular flow.

Figure 1 Overview of laboratory measurements of Cappa et al. (2019). The fault slip and dilatancy recorded over the full creep stage of the experiment are shown in (a) and (c), respectively, along with the fluid pressure for reference. The final stage of the experiment (grey-shaded area of a and c) is enlarged in (b) and (d).

Figure 2 Lower triangle (blue graphs and those below them): corner plot of the posterior distributions of the inverted parameters L, H, ϕ₀, and ϕ_c (marginalised over the nuisance parameter ν²) for the laboratory injection experiment. The main diagonal panels show the posterior probability density distribution of each parameter, whereas the off-diagonal panels show the co-variance of posterior samples. The black dashed lines in the distribution plots (main diagonal) indicate the median value. Upper triangle (green graphs and the one above them): corner plot of the posterior distributions of A and B (see main text).

We first fit Eq. (6) to the measured dilatancy as a function of slip. Since the data are sampled uniformly in time but not in slip (as the sample deformation is accelerating), we interpolate the slip data to assign uniform weight to each measurement during the inversion. The bounds on the prior distribution are given by $\mathrm{10}<L<\mathrm{100}$ µm, $\mathrm{0.1}<H<\mathrm{1}$, $\mathrm{0.1}<{\mathit{\varphi }}_{\mathrm{0}}<{\mathit{\varphi }}_{\mathrm{c}}$, and $\mathrm{0.2}<{\mathit{\varphi }}_{\mathrm{c}}<\mathrm{0.4}$. The resulting posterior distributions are presented in a corner plot (Fig. 2), showing appreciable trade-offs between L and H and between ϕ₀ and ϕ_c. Nonetheless, the parameters L and H are reasonably well resolved as $L=\mathrm{64.3}\pm \mathrm{8.9}$ µm and $H=\mathrm{0.45}\pm \mathrm{0.06}$ (median ± 1 SD). And while ϕ₀ and ϕ_c trade off almost perfectly and hence span a near-uniform distribution over the permitted parameter range, their difference is well resolved as ${\mathit{\varphi }}_{\mathrm{c}}-{\mathit{\varphi }}_{\mathrm{0}}=\mathrm{0.09}\pm \mathrm{0.01}$. These parameter values are perfectly consistent with previous studies (e.g. Chen and Spiers, 2016; van den Ende et al., 2018). Although the inferred layer thickness L is much less than the total thickness of the sample (initially around 5 mm), one should keep in mind that deformation localises in a much narrower zone, so that the effective thickness of the actively deforming region of the gouge is much less than the total sample thickness. In similar experiments conducted by Scuderi et al. (2017), the localised region was observed to have a thickness of 10–20 µm, which was inevitably affected by post-experiment compaction. Moreover, the experiment was performed in a double-direct shear configuration, so that the total thickness inferred here represents the thickness of localised regions on both sides of the central forcing block. In other words, the total layer thickness inferred from the post-mortem microstructural observations would be at least 20–40 µm. Hence, our inferred estimate of 64 µm seems appropriate for an actively deforming localised gouge layer. Niemeijer and Spiers (2007) derived a theoretical lower bound on $H\ge \sqrt{\mathrm{3}}/\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}\approx \mathrm{0.96}$ for a dense 2D packing of hexagonal grains. Since the third dimension plays an important role in strain accommodation within granular media (Frye and Marone, 2002; Hazzard and Mair, 2003), we expect this lower bound to be lower in a 3D system with more degrees of freedom. Upon inspection of Eq. (6), we can formulate the mapping between layer thickness and slip as $\mathrm{\Delta }L=A\left[\mathrm{1}-\exp \left(-\mathrm{2}B\mathit{\delta }\right)\right]$ and infer $A=L\left({\mathit{\varphi }}_{\mathrm{c}}-{\mathit{\varphi }}_{\mathrm{0}}\right)/\left(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{c}}\right)$ and $B=H\left(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{c}}\right)/L$ as lumped parameters (upper triangle of Fig. 2). Since A and B are the only parameters directly constrained by the data, the original four parameters depend on them and show strong trade-offs.

Figure 3 Panel (b) of the posterior distributions of the inverted parameters $\stackrel{\mathrm{̃}}{a}$ and $\stackrel{\mathrm{̃}}{\mathit{\mu }}$ (marginalised over the nuisance parameter ν²) for the laboratory injection experiment. Panels (a) and (c) show the posterior probability density distribution of each parameter, whereas (b) shows the co-variance of posterior samples. The black dashed lines in the distribution plots (a, c) indicate the median value.

We continue by fitting the (logarithm of) measured slip rate based on Eq. (3), using the parameter values inferred in the previous step to compute the time evolution of tan ψ. Without loss of generality, we define ${\dot{\mathit{\gamma }}}_{\text{gr}}^{*}=\mathrm{1}$ µm s⁻¹ L⁻¹, so that ${\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{*}$ represents the grain boundary friction coefficient at a slip rate of $V=\mathrm{1}\mathrm{\mu }\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. Since ϕ₀ and ϕ_c individually are ambiguous, we take ϕ₀=0.25 and increment this value by the inverted ϕ_c−ϕ₀ to obtain ϕ_c=0.34. The slip rate parameters are extremely well resolved (see Fig. 3), and found to be $\stackrel{\mathrm{̃}}{a}=\left(\mathrm{10.26}\pm \mathrm{0.15}\right)\times {\mathrm{10}}^{-\mathrm{3}}$ and ${\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{*}=\mathrm{0.6852}\pm \mathrm{0.00028}$, with minimal trade-off between the two parameters. With these parameters, the fit to the slip rate data is excellent (Fig. 4b).

Figure 4 Results of the inversion procedure. (a) Dilatancy versus slip, with the inversion curve given by Eq. (6); (b) slip rate versus time, with the inversion curve given by Eq. (1a); (c, d) forward model results of slip and dilatancy for a (narrow) range of initial porosity ϕ₀ (as indicated in the legend). The reference value of ϕ₀ is obtained from the inversion of Eq. (6).

Finally, for verification, we numerically solve the forward model given by Eqs. (1a) and (1b) with the parameters obtained in the inversion procedure (Fig. 4c and d). While we obtain an excellent fit with the observed time evolution of slip and dilatancy, we also find that the forward model is extremely sensitive to the initial condition ϕ₀. While the overall features of the simulated sample response are similar, the exponential sensitivity to porosity leads to critical behaviour and strong variations in the timing of the sample failure. Since the rate of increase in porosity is proportional to the shear strain rate, which in turn is an exponential function of porosity (refer to Eqs. 1b and 3), the positive feedback loop leads to an extremely rapidly diverging state. This is highlighted in Fig. 4c and d, where we vary the initial porosity between 0.2491 and 0.2514 (−1.5 % and +1 % around the reference value of 0.25). The initial condition that gives the best match in terms of the onset of accelerated slip is close to the initially chosen value of ϕ₀=0.25, although we assign no significance to tiny deviations in the initial porosity. In a laboratory setting, the sensitivity of the modelled slip rate falls well within the measurement resolution of the sample porosity (typically of the order of several percent of units of porosity), so verification of this sensitivity would be challenging.

3.2 Field experiment

Encouraged by the results of the proposed inversion method for the laboratory experiment, we continue to apply the same procedure to the field injection test of Cappa et al. (2019). Like in the laboratory experiment, the in situ pressurisation of a tectonic fault triggered accelerating slip and, associated with it, fault opening (dilatancy) – see Fig. 5a. While the acceleration of shear and normal displacement on the fault was more gradual than in the laboratory experiment, a phase of rapidly accelerating slip at t>800 s can be clearly seen. The amount of dilatancy measured as a function of slip (Fig. 5b) was proportionally more than in the laboratory experiment by at least 1 order of magnitude, so we expect a priori that the frictional parameters inferred from the laboratory cannot immediately describe the behaviour of the fault in situ. Indeed, when we perform the inversion of the dilatancy-slip data from the field experiment, we find median values of $L=\mathrm{4.9}\pm \mathrm{1.0}\mathrm{mm}$, $H=\mathrm{13.0}\pm \mathrm{2.5}$, and ${\mathit{\varphi }}_{\mathrm{c}}-{\mathit{\varphi }}_{\mathrm{0}}=\mathrm{0.082}\pm \mathrm{0.016}$ (Fig. 6).

Figure 5 (a) Measurements of fault shear and normal displacement and fluid pressure during the field injection test of Cappa et al. (2019); (b) inversion of the dilatancy measured during the field injection experiment. To produce a reasonable fit to the data, an unrealistically high value of the dilatancy parameter H=6.4 was required.

Figure 6 Corner plot of the posterior distributions of the inverted parameters L, H, ϕ₀, and ϕ_c (marginalised over the nuisance parameter ν²) for the field injection experiment. The main diagonal panels show the posterior probability density distribution of each parameter, whereas the off-diagonal panels show the co-variance of posterior samples.

While these other values seem entirely reasonable, the inferred value of H is well above the estimated upper bound of $\sqrt{\mathrm{3}}/\mathrm{2}{\mathit{\varphi }}_{\mathrm{c}}\approx \mathrm{2.9}$. This suggests that the CNS model is unable to explain the relationship between fault slip and fault opening in this experiment. In the CNS model, dilatancy is envisioned to originate from grain sliding and rolling, neighbour swapping, and “jostling”, which requires a volume increase of the gouge to accommodate. However, the model fault itself is mathematically planar, and so no dilatancy occurs due to geometric constraints. In the case of a macroscopically non-planar fault geometry (as is inevitable for tectonic faults; Candela et al., 2012), additional dilatancy (with associated permeability changes) at the onset of slip is necessary. Moreover, poro-elastic effects (elastic fault opening) due to fluid pressure changes are not considered here. The inability of the CNS model to describe the fault opening with a reasonable choice of parameters is therefore not a shortcoming of the CNS model (which describes the mechanics of a small representative volume element), but is rather due to an incomplete coupling with processes that transcend the scale envisioned by the CNS model. For simple, spatially uniform relationships between geometric fault opening and fault slip, this first-order contribution to the fault dilatation may be incorporated into Eq. (1b). However, for more realistic (i.e. spatially heterogeneous) fault opening, a multi-scale numerical extension of the adopted model is required.

Since the CNS model fault strength (and therefore the fault slip rate) is directly controlled by the dilatancy parameter H, it is unwarranted to attempt to infer $\stackrel{\mathrm{̃}}{a}$ and ${\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{*}$ based on the parameters inferred from the dilatancy. While this may seem like a severe limitation of the CNS model, it actually serves as an important indication of the applicability of the model, and the validity of its parameters, when attempting to extrapolate to nature. Moreover, the basic mechanics of the CNS model are still expected to govern the strength and slip rate of the fault, even though part of the model predictions (the dilatancy) cannot be constrained by independent measurements. By numerically solving the forward model, the fault slip as a function of time and fluid pressure may be reproduced within a reasonable range of parameter values, for which the predicted fault opening would likely be much less than measured by Cappa et al. (2019).

4 Discussion

4.1 Comparison with rate-and-state friction

Traditionally, laboratory experiments are interpreted within the framework of rate-and-state friction (RSF), commonly presented as (Dieterich, 1979; Ruina, 1983)

$$\begin{array}{}\text{(7a)}& {\displaystyle }\mathit{\mu }\left(V,\mathit{\theta }\right)& {\displaystyle }={\mathit{\mu }}^{*}+a\ln \left({\displaystyle \frac{V}{V^{*}}}\right)+b\ln \left({\displaystyle \frac{V^{*}\mathit{\theta }}{D_{\mathrm{c}}}}\right),\text{(7b)}& {\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}t}}& {\displaystyle }=\left\{\begin{array}{ll}\mathrm{1}-\frac{V\mathit{\theta }}{D_{\mathrm{c}}},& \text{“ageing law”}\\ -\frac{V\mathit{\theta }}{D_{\mathrm{c}}}\ln \left(\frac{V\mathit{\theta }}{D_{\mathrm{c}}}\right),& \text{“slip law”}\end{array}\right.,\end{array}$$

where μ(V,θ) is the macroscopic friction coefficient at slip rate V and “state” θ, μ^* is a reference friction coefficient at slip rate V^*, and a, b, and D_c are empirical constants. As has been shown by Chen et al. (2017), the CNS model is asymptotically identical to RSF for small departures from steady state, for which the CNS equivalents of the RSF parameters a, b, and D_c can be treated as constants. For large departures from steady state, the behaviours predicted by CNS and RSF diverge, as the aforementioned parameters can no longer be considered to be constant (van den Ende et al., 2018). Nonetheless, because of their similarity, the limitations of the CNS model also apply to rate-and-state friction. One advantage of using the CNS microphysical model over traditional RSF, is that the governing parameters have a more physically meaningful interpretation. Even though numerous studies have attempted to elucidate the physical origin of RSF (Aharonov and Scholz, 2018; Brechet and Estrin, 1994; Ikari et al., 2016; Putelat et al., 2011), in practice these theoretical constraints are not considered. Instead, it is more convenient to constrain the RSF parameters empirically through laboratory velocity-step experiments (Blanpied et al., 1998; Carpenter et al., 2016; Chester, 1994; Hunfeld et al., 2017; Reinen and Weeks, 1993). With these laboratory measurements of the RSF parameters, fault slip observed during decametric-scale fluid injection tests can be accurately modelled (Cappa et al., 2019), although the same behaviour can be obtained for a wide range of parameter values: in the study of Cappa et al. (2019) a similar fit to the data was obtained for velocity-weakening ($\left(a-b\right)<\mathrm{0}$) and velocity-strengthening ($\left(a-b\right)>\mathrm{0}$) friction, even though seismic slip can only be produced in the former case of velocity-weakening friction. Hence, more observational constraints are required to distinguish between the different types of behaviour.

Aside from the fault-parallel slip, the fault opening potentially provides a second prominent constraint. In the classical RSF framework, volumetric deformation is not explicitly accounted for. Traditionally, the state parameter θ has been interpreted as encoding the average lifetime of asperity contacts (at steady state) or the relative area of asperity contacts (Dieterich, 1994; Scholz, 2019), both of which do not entail volumetric deformation of the fault gouge. Empirical relations between the state parameter θ and porosity have also been proposed (Segall and Rice, 1995; Sleep, 2005) and used in hydro-mechanical modelling (Jeanne et al., 2018), but these relations are typically not employed as additional constraints of the RSF constitutive parameters. Moreover, relations between the steady-state coefficient of friction (and its velocity dependence) have been established based on energy balance considerations (Beeler et al., 1996; Marone et al., 1990). Since these relations pertain to the steady-state coefficient of friction, they do not apply to non-steady-state conditions (for which $\mathrm{d}V/\mathrm{d}t\ne \mathrm{0}$ and $\mathrm{d}\mathit{\varphi }/\mathrm{d}t\ne \mathrm{0}$) and do not offer additional insight on the relationship between θ and ϕ. On the other hand, volumetric deformation is an integral part of the CNS model, hence allowing (and requiring) us to incorporate these measurements to arrive at a better constrained set of parameters.

4.2 Relationships between experiments and nature

While the CNS microphysical parameters can be directly estimated from laboratory experiments, their incorporation into numerical models of tectonic faults may be subject to moderation based on geological or physical considerations. In laboratory experiments conducted at room ambient conditions and comparatively high deformation rates (of the order of micrometres per second up to millimetres per second), the gouge porosity remains close to the critical-state porosity. Likewise, in the laboratory experiment of Cappa et al. (2019), the initial porosity was estimated to be less than 0.1 units of porosity below the critical-state porosity. Given longer timescales and higher temperatures, compaction induced by one or more time-dependent creep mechanisms (such as pressure solution creep or subcritical crack growth) would gradually reduce the porosity of the gouge, thereby increasing its strength and critical fluid pressure at which the fault slip rates become appreciable. In numerical simulations of fault slip, the initial state of a tectonic fault is likely not the same as for the laboratory fault. Fortunately, this initial state could be estimated from microstructural analyses of drill cores. Moreover, the choice of the initial state of the fault does not affect any of the other frictional parameters of the CNS model. This is in contrast to rate-and-state friction, where the initial value of the state parameter (θ at t=0) should also affect the magnitude of b, which has been found to increase with decreasing porosity (or increasing θ; Chen et al., 2015, 2017).

The property that b (or more precisely: b∕D_c) is sensitive to the gouge porosity can also be derived from stability analysis of the CNS model. Consider the general criterion for unstable slip of a spring block:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}\mathit{\tau }}{\mathrm{d}\mathit{\delta }}}={\displaystyle \frac{\mathrm{1}}{V}}\left({\displaystyle \frac{\partial \mathit{\tau }}{\partial \mathit{\varphi }}}{\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}t}}+{\displaystyle \frac{\partial \mathit{\tau }}{\partial V}}{\displaystyle \frac{\mathrm{d}V}{\mathrm{d}t}}\right)\le -K,\end{array}$$

where K is the shear stiffness of the fault. For an instantaneous step change in velocity, $\mathrm{d}V/\mathrm{d}t=\mathrm{0}$ for t>0. Assuming that unstable slip is governed by the onset of granular flow, the shear strength is given by the CNS model as (Chen and Spiers, 2016)

$$\begin{array}{}\text{(9)}& \mathit{\tau }={\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathit{\mu }}+\tan \mathit{\psi }}{\mathrm{1}-\stackrel{\mathrm{̃}}{\mathit{\mu }}\tan \mathit{\psi }}}{\mathit{\sigma }}_{\mathrm{e}}.\end{array}$$

Hence, using Eq. (1b), the stability criterion can be expressed in terms of microstructural quantities as (van den Ende et al., 2018)

$$\begin{array}{}\text{(10)}& K\le \mathrm{2}H\left(\mathrm{1}-\mathit{\varphi }\right)\tan \mathit{\psi }{\displaystyle \frac{\mathrm{1}+{\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{\mathrm{2}}}{{\left(\mathrm{1}-\stackrel{\mathrm{̃}}{\mathit{\mu }}\tan \mathit{\psi }\right)}^{\mathrm{2}}}}{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{e}}}{L}}.\end{array}$$

In the vicinity of steady state, the above statement should be identical to the stability criterion derived from rate-and-state friction, i.e. (Rubin and Ampuero, 2005)

$$\begin{array}{}\text{(11)}& K\le K_{\mathrm{b}}={\displaystyle \frac{b{\mathit{\sigma }}_{\mathrm{e}}}{D_{\mathrm{c}}}}.\end{array}$$

Here, K_b is a critical stiffness value that facilitates acceleration of slip (seismic or aseismic). From the comparison of the two inequalities, it can be concluded that the K_b therefore must increase with decreasing porosity. This was also observed in the discrete element model simulations of van den Ende and Niemeijer (2018), which were conducted completely independently of the assumptions and limitations of the CNS model. We note that the comparison between Eqs. (10) and (11) only holds in the vicinity of steady state. Nonetheless, Eq. (10) can be used to describe the stability of fault slip far from steady state, circumventing the issue of the velocity and state dependence of a, b, and D_c (as observed by Cappa et al., 2019; den Hartog and Spiers, 2013; Reinen et al., 1992; Takahashi et al., 2017, and many others). For a more detailed analysis of the frictional stability of a model fault governed by the CNS model, we refer to Chen and Niemeijer (2017).

Combining now the observations made in Sect. 3 with the discussion above, we propose that the seismogenic potential of faults subjected to fluid pressure perturbations is best described in terms of the dilatant behaviour of the fault and its initial degree of compaction. One can infer the microphysical parameters H, ϕ_c, $\stackrel{\mathrm{̃}}{a}$, and ${\stackrel{\mathrm{̃}}{\mathit{\mu }}}^{*}$ from laboratory experiments, and assume reasonable in situ values of ϕ and L for the tectonic fault to simulate its response to a changing stress field (fluid pressure). If permitted by the numerical method, fault non-planarity, permeability changes, and elastic moduli reduction may be introduced to add further complexity, as anticipated based on the results of Sect. 3.2. In this way, the evident pressure and velocity dependence of the rate-and-state friction parameters can directly be accounted for in a self-consistent and transparent manner and the model outcomes interpreted in terms of physical observables.

5 Conclusions

In this work, we analysed the fluid injection experiments conducted by Cappa et al. (2019) in the laboratory and in situ, in terms of the Chen–Niemeijer–Spiers (CNS) microphysical model. We proposed a Bayesian inversion approach to extract the governing parameters without the need for numerically solving the forward problem, while elucidating the uncertainties and trade-offs between the model parameters. We showed that while the localised gouge layer thickness L and the dilatancy parameter H can be well resolved, the initial and critical-state porosities trade off perfectly, so that only their difference ϕ_c−ϕ₀ can be resolved in the experiments. When numerically solving the forward model with the inferred parameter values, we obtained almost perfect agreement with the measurements, indicating that the CNS model accurately describes fault deformation in response to a fluid pressure perturbation. When the same inversion approach was applied to a decametric-scale field injection experiment, we found that the inferred parameters fell outside of the feasible range of values, highlighting the relevance of other mechanisms, such as fault structure and poro-elastic effects, in this scenario.

The excellent agreement between the CNS model and the laboratory data allows us to interpret the dynamics of the fault in terms of volumetric deformation (porosity changes). By doing so, we circumvent the velocity dependence of the rate-and-state friction parameters a, b, and D_c, which increases the predictive power of numerical models of natural faults. Adopting the CNS model expedites the extrapolation of laboratory results to nature and permits better assessment of the applicability of the model and accuracy of the parameter values.