TY - JOUR
AU - Gesret,, A
AB - SUMMARY In situ observations of fluid induced fault slip reactivation, as well as the analysis of induced seismicity have demonstrated the complexity of fluid–fault interactions under geological conditions. If fluid flow commonly reactivates faults in the form of aseismic slip or earthquakes, the resulting shear deformation causes strong modifications of the hydraulic properties. In this context, the relationship between slip front and fluid front on deep faults remains not fully understood. In this study, we investigate shear induced fluid flow and hydraulic diffusivity enhancement during fracture shearing in the laboratory. We use a series of injection reactivation tests, conducted under triaxial conditions, at different confining pressures (30, 60 and 95 MPa). The evolution of the fluid pressure along the saw-cut Andesite rock sample was monitored by two pressure sensors, at two opposite locations of the experimental fault. We estimate the history of the effective hydraulic diffusivity (and its associated uncertainties) governing the experimental fault, using the pressure history at two points on the fault. For this, we develop a deterministic and a probabilistic inversion procedure, which is able to reproduce the experimental data for a wide time range of the different experiments. In this study, the hydraulic diffusivity increases by one order of magnitude through the injection experiment. Hydraulic diffusivity changes are mainly governed by the reduction of the effective normal stress acting on the fault plane, with a second-order effect of the shear slip. Fracture and flow, Inverse theory, Induced seismicity 1 INTRODUCTION Understanding how the permeability of a fault evolves during fluid injection activities and through the fault reactivation process is of great interest, for reservoir engineering, enhanced geothermal systems, as well as for hydraulic fracturing operations. It can also help better understand the spatio-temporal distribution of induced seismic sequences. However, the interactions between fluids and faults/reservoirs evolve in time and space, as the confining pressure, effective stress and shear slip can affect the hydromechanical properties of the fault. On the one hand, variation in fault permeability has been observed following changes in effective stress during laboratory experiments (Zoback & Byerlee 1975; McKee et al. 1988; Ghabezloo et al. 2009; Rutter & Mecklenburgh 2018) as well as during in situ permeability measurements (Fisher & Zwart 1996). On the other hand, seismic events and slip accumulation can also affect the fault permeability (Zhang & Tullis 1998; Baghbanan & Jing 2008), as permeability enhancement was observed similarly during laboratory (Chen et al. 2000a; Gutierrez et al. 2000; Wu et al. 2017; Im et al. 2018) and in situ injection experiments (Guglielmi et al. 2015a, b; Duboeuf et al. 2017; Bhattacharya & Viesca 2019). Current permeability measurement methods consist mainly of experimental approaches, in which the permeability is determined via fluid flow analysis, using Darcy’s law for instance (Darcy 1857), either on experimental samples in the laboratory (e.g., Zoback & Byerlee 1975; Zhang & Tullis 1998; Ghabezloo et al. 2009) or in situ along a plate-boundary fault (e.g. Fisher & Zwart 1996). Although quite effective, such classical methods lie on the assumption that the permeability is constant throughout one measurement test. Thus, they do not allow for the characterization of the permeability evolution resulting from slip or changes in effective stress, in the context of a single laboratory injection experiment per se. Beyond permeability enhancement, the relation between the pressure diffusive front and the fault reactivation process is the focus of several research studies. The migration of seismic events has been suspected to be driven by the diffusion of pore pressure away from the injection wells (Shapiro et al. 1997). This is proposed for instance for the induced seismic sequences in Denver, Colorado (Healy et al. 1968; Hermann et al. 1981; Hsieh & Bredehoeft 1981), Soultz-sous-Forêt (Shapiro et al. 2002) and Ohio (Kim 2013). This mechanism assumes that the reactivation front tracks a particular pressure front. In addition, a recent induced seismicity triggering mechanism has been proposed by De Barros et al. (2018), where authors argue that aseismic deformation could trigger seismic ruptures through stress transfer beyond the pressurized region. A number of numerical and theoretical studies support this theory and predict that the aseismic deformation front can in some cases outpace fluid diffusion (Garagash & Germanovish 2012; Galis et al. 2017; Bhattacharya & Viesca 2019; Dublanchet 2019); however no direct observations have been made so far, as it is difficult to trace the diffusion front for real cases of induced seismicity. Understanding the relationship between slip (reactivation) front and fluid front requires a better understanding of what controls hydraulic diffusivity in this context. We thus propose a numerical method, in the context of deterministic and probabilistic inversion approaches, that allows to estimate the temporal evolution of the hydraulic diffusivity of an experimental fault throughout an injection test, using the pressure history at two points on the fault. The deterministic approach is a gradient-based approach, used to determine the optimal model in a least-squares sense, and the probabilistic one is used to evaluate the associated uncertainties. In this study, we investigate shear induced fluid flow and permeability enhancement during fracture shearing. We used a series of laboratory injection tests on saw cut Andesite rock sample, under triaxial conditions, in which water was injected under constant rate into the experimental fault. The evolution of the fluid pressure along the saw-cut Andesite rock sample was monitored by two pressure sensors, at two opposite locations of the experimental fault, while the evolution of the strain was observed using strain gages along fault strike. We then developed and implemented a gradient-based approach (deterministic approach), with the use of the adjoint state technique (Plessix 2006), and a Markov Chain Monte Carlo algorithm (probabilistic approach, Robert & Casilla 2004) in an inversion framework to the experimental data. The objective is to estimate the hydraulic diffusivity enhancement and interpret it with respect to the accumulated shear slip and the reduction in effective stress. In the following, we present the experimental setup and data, then introduce the numerical methods and their application to the data, to finally discuss our numerical findings. 2 EXPERIMENTAL DATA 2.1 Experimental Setup Triaxial shear experiments were conducted on an Andesite cylindrical rock sample, which presents a young modulus E = 64 GPa and a porosity ϕ = 2.0 per cent (Li et al. 2019). The cylinder has a length of H = 8.8 cm and a radius of R = 2.0 cm (values in Table 1). A sketch of the front view of the experimental setup is presented in Fig. 1. A single saw cut fracture represents the experimental fault. It forms a 30° angle with respect to the vertical axis of the rock sample, and has an elliptical shape of 8.0 cm length along strike and 4.0 cm wide (values in Table 1). Two vertical boreholes were drilled up to the fault surface, starting from the bottom and top of the rock sample, respectively, at opposing edges of the elliptical fault plane (Fig. 1). The boreholes were drilled 5.0 mm away from the edge of the rock sample, and have a 4.0 mm diameter. The bottom borehole was connected to the pump A, where liquid water was injected, and served as the injection borehole; while the top one was connected to the pump B and was sealed during the injection experiment, and served as an observation borehole. The experimental sample was equipped by 8 strain gages that were uniformly distributed along the fault plane, as shown in Fig. 1. Figure 1. Open in new tabDownload slide Sketch of the experimental setup: H is the cylinder length and R its radius, L is the fault length, α its angle with respect the the long axis of the cylinder, σ is the axial loading and Pc is the confining pressure. Figure 1. Open in new tabDownload slide Sketch of the experimental setup: H is the cylinder length and R its radius, L is the fault length, α its angle with respect the the long axis of the cylinder, σ is the axial loading and Pc is the confining pressure. Table 1. Laboratory injection experiment: list of geometrical and mechanical parameters. Parameter . Symbol . Value . Unit . Cylinder length H 8.8 cm Cylinder radius R 2 cm Fault length L 8 cm Fault width W 4 cm Fault angle w.r.t σ α 30 (°) Young modulus E 64 GPa Porosity ϕ 2 per cent Parameter . Symbol . Value . Unit . Cylinder length H 8.8 cm Cylinder radius R 2 cm Fault length L 8 cm Fault width W 4 cm Fault angle w.r.t σ α 30 (°) Young modulus E 64 GPa Porosity ϕ 2 per cent Open in new tab Table 1. Laboratory injection experiment: list of geometrical and mechanical parameters. Parameter . Symbol . Value . Unit . Cylinder length H 8.8 cm Cylinder radius R 2 cm Fault length L 8 cm Fault width W 4 cm Fault angle w.r.t σ α 30 (°) Young modulus E 64 GPa Porosity ϕ 2 per cent Parameter . Symbol . Value . Unit . Cylinder length H 8.8 cm Cylinder radius R 2 cm Fault length L 8 cm Fault width W 4 cm Fault angle w.r.t σ α 30 (°) Young modulus E 64 GPa Porosity ϕ 2 per cent Open in new tab The rock sample was subjected to an axial loading σ and radial one Pc (or confining pressure). The normal σn and shear stress τn (average values on the fault plane) can be estimated by projecting the triaxial stress state onto the fault plane, as follows: $$\begin{eqnarray*} \sigma _\mathrm{n} = \left( \frac{\sigma +P_\mathrm{c}}{2} \right) - \left( \frac{\sigma - P_\mathrm{c}}{2}\right) \mathrm{cos}(2 \alpha ), \end{eqnarray*}$$(1) $$\begin{eqnarray*} \tau _\mathrm{n} = \left( \frac{\sigma - P_\mathrm{c}}{2}\right) \mathrm{sin}\left( 2 \alpha \right). \end{eqnarray*}$$(2) The two boreholes allow us to get a continuous measure of the pressure at two locations along the fault plane (p|$_\mathrm{inj}^E$| and p|$_\mathrm{obs}^E$|⁠) throughout the injection experiment. Average fault slip can also be computed by projecting the axial displacement onto the fault plane, and strain deformations were recorded using the strain gages. The fluid flow parameters (pressure, volume and injection rate) were recorded at only 1 Hz, estimated to be sufficient to track the diffusion process. 2.2 Experimental protocol We carried out experiments at different confining pressures (Pc = 30, 60 and 95 MPa). At the start of the injection experiments, the pore pressure was set to 10 MPa uniformly along the fault plane. Then a loading phase was conducted: the shear stress was increased to ∼90 per cent of the peak shear stress (τ0 = 0.9τp). For each injection experiment, the peak shear stress was determined by conducting a prior axial loading test. At the end of the loading phase, liquid water was injected throughout the injection borehole, under a constant injection pressure rate of ∼5 MPa min–1, until reaching a pre-set value of maximum pressure (lower than the confining pressure). Injection continued until the pressure equilibrium is reached along the fault plane. As liquid water was injected at ambient temperature, we consider the injection isothermal and neglect any thermo-elastic behaviour. Values of the different injection experiment parameters are listed in Table 2. Table 2. Laboratory injection experiment: list of experimental parameters. Parameter . Symbol . Exp. 1 . Exp. 2 . Exp. 3 . Unit . Injection rate 5 5 5 MPa min–1 Initial pore pressure pinit 10 10 10 MPa Confining pressure Pc 30 60 95 MPa Maximum pressure pmax ≈29 ≈59 ≈89 MPa Experiment duration 1063 1877 5548 s Injection start time 143 147 670 s Parameter . Symbol . Exp. 1 . Exp. 2 . Exp. 3 . Unit . Injection rate 5 5 5 MPa min–1 Initial pore pressure pinit 10 10 10 MPa Confining pressure Pc 30 60 95 MPa Maximum pressure pmax ≈29 ≈59 ≈89 MPa Experiment duration 1063 1877 5548 s Injection start time 143 147 670 s Open in new tab Table 2. Laboratory injection experiment: list of experimental parameters. Parameter . Symbol . Exp. 1 . Exp. 2 . Exp. 3 . Unit . Injection rate 5 5 5 MPa min–1 Initial pore pressure pinit 10 10 10 MPa Confining pressure Pc 30 60 95 MPa Maximum pressure pmax ≈29 ≈59 ≈89 MPa Experiment duration 1063 1877 5548 s Injection start time 143 147 670 s Parameter . Symbol . Exp. 1 . Exp. 2 . Exp. 3 . Unit . Injection rate 5 5 5 MPa min–1 Initial pore pressure pinit 10 10 10 MPa Confining pressure Pc 30 60 95 MPa Maximum pressure pmax ≈29 ≈59 ≈89 MPa Experiment duration 1063 1877 5548 s Injection start time 143 147 670 s Open in new tab 2.3 Experimental results Fig. 2 displays the measurements of the pore pressure (at injection and observation boreholes) and the average slip throughout the three injection experiments. We should note that the instantaneous decrease then increase in pressure at 2500 s for the test at 95 MPa confining pressure (Fig. 2c) is caused by a slight imperfection in the experimental protocol: pump A was emptied and replaced. Figure 2. Open in new tabDownload slide Experimental gross results: time-series of pore pressure (left-hand axis) and average slip along strike (right-hand axis) for the different injection experiments at different confining pressure values Pc: (a) 30 MPa, (b) 60 MPa, here the small subplot is a zoom around the initial time of the injection and (c) 95 MPa. In all subplots, the vertical dashed line represents the limit between the loading phase and the injection phase. Points A–I represent examples of slip events. Figure 2. Open in new tabDownload slide Experimental gross results: time-series of pore pressure (left-hand axis) and average slip along strike (right-hand axis) for the different injection experiments at different confining pressure values Pc: (a) 30 MPa, (b) 60 MPa, here the small subplot is a zoom around the initial time of the injection and (c) 95 MPa. In all subplots, the vertical dashed line represents the limit between the loading phase and the injection phase. Points A–I represent examples of slip events. Upon the start of injection, instantaneous increases in pressure at the observation borehole are observed, especially for the two injection tests at 30 and 60 MPa confining pressure. These increases can not be due to pure diffusion effects and may be due to some direct poro-elastic effects. As in this study we only model the fluid diffusion process, we will not interpret the instantaneous response for the rest of the study. From Fig. 2, we observe that fluid injection and pressure diffusion reactivate the experimental fault. Various slip events are observed during the injection phase for the different tests. We selected in particular a few examples (events A–I) to illustrate the existence of different types of slip events: for instance events A, B, D, E and G are accompanied by a large instantaneous slip along the fault plane, while the rest does not share this feature. These slip events are characterized by relatively higher average slip rates with respect to the other events. Note that the average slip rate is estimated as the time derivative of the average slip measured along the fault. In particular, the average slip rate reaches 0.1 mm s–1 for event A, 0.69 mm s–1 for event D and 7.48 mm s–1 for event G. The slip events are associated with instantaneous pressure drops in the injection borehole, and strong increases of pressure in the observation borehole. These changes are particularly significant during events A and B, D and E and G and H at 30, 60 and 95 MPa confining pressure, respectively. These instantaneous pressure changes could be interpreted by an enhanced pressure diffusion, that is an increase in the hydraulic diffusivity. To summarize, the results show that the experimental fault can be reactivated due to fluid injection through different slip events, and may be accompanied by an enhancement of hydraulic diffusivity. The latter is the main focus of this study. In order to better understand the role of the different parameters at play, we aim to characterize the diffusivity history throughout the injection experiment, as it is suspected to vary with the effective stress reduction, as well as following the different slip events recorded on the fault. We perform numerical inversion on the experimental data in a deterministic and probabilistic framework. In the next two sections, we present the development of the numerical method and its application to the experimental data. 3 METHODOLOGY The understanding of shear induced diffusivity enhancement requires to invert the spatio-temporal diffusivity history D(x,y,t) along the fault throughout the injection experiment, from the observed pressure p|$_\mathrm{obs}^E$|(t) (we note that the superscripts E and N here stand for experimental and numerical, respectively). The direct problem involves the pressure history p(x,y,t) and is given by: $$\begin{eqnarray*} \frac{\partial {p}}{\partial {t}} - {D}({t}) \Delta {p} = 0, \end{eqnarray*}$$(3) where Δ is the Laplacian operator, p(x, y, t = 0) = 10 MPa, |${p(x=x_{inj},y=y_{inj},t)} = {p_{inj}^E(t)}$| and ∂p/∂n = 0 (n being the normal direction) at the fault boundaries representing no fluid flow. As we have only one observation position, it would be highly under-determined to inverse for the full spatio-temporal history of the diffusivity. Thus, we choose to inverse for an effective diffusivity D(t) that only depends on time in order to reduce the number of parameters to be determined. For that, we first perform a deterministic inversion in order to find the diffusivity history that minimizes the distance to the observed pressure. In a second step, we conduct a probabilistic inversion to determine the uncertainty of our diffusivity history. 3.1 Deterministic inversion The deterministic inversion is performed using the adjoint state method (Plessix 2006). Let λ(x, y, t) be the adjoint variable, the adjoint problem is given by: $$\begin{eqnarray*} \frac{\partial \lambda }{\partial {t}} + {D}({t})~\Delta \lambda = \delta ({x - x_{obs}}) \delta ({y - y_{obs}}) \left( {p_{obs}^N} - {p_{obs}^E} \right), \end{eqnarray*}$$(4) where λ(x, y, t = Tmax) = 0, and ∂λ/∂n = 0 at the fault boundaries). In this case, the objective function J reads: $$\begin{eqnarray*} J[D,\lambda ,{p}] &=& \frac{1}{2} \iiint \left( \delta \left( {x} - {x}_\mathrm{obs} \right) \delta \left( {y} - {y}_\mathrm{obs} \right) \left( {p}_\mathrm{obs}^N - {p}_\mathrm{obs}^E \right)^2 \right) {{\rm d}x~{\rm d}y~{\rm d}t} \\ && - \, \iiint \lambda ({x}, {y}, {t}) . \left[ \frac{\partial {p}}{\partial {t}} - {D} \Delta {p} \right] {{\rm d}x~{\rm d}y~{\rm d}t}. \end{eqnarray*}$$(5) The gradient ∂J/∂D(t) is the essential element and represents the correction or the update of the model: $$\begin{eqnarray*} \frac{\text{d} {J}}{\text{d} {D} ({t}_0)} = \int \left( \delta \left( {t} - {t}_0 \right) \left[ \frac{ \partial {p}_\text{obs}^\text{N}({t})}{\partial {D}({t}_0)} \right].\left[ {p}_\text{obs}^\text{N}({t}) - {p}_\text{obs}^\text{E}({t}) \right] \right) {{\rm d}t}. \end{eqnarray*}$$(6) We then make use of a basic gradient descent algorithm to optimize the gradient given by: $$\begin{eqnarray*} \frac{{{\rm d}J}}{{{\rm d}D}} = \iint \lambda ~\Delta {p}~{{\rm d}x~{\rm d}y}. \end{eqnarray*}$$(7) A summary of the adjoint state method to the 2-D pressure diffusion problem in the context of this study is presented in Appendix A. This process is iterated until the reduction of the objective function [J(n) – J(n+1)]/J(n) becomes smaller than 10−3. The two differential eqs (3 and 4) are solved using an explicit time scheme of the finite volume method. We consider a very simple discretization grid to the fault plane with equal x and y spatial discretization: Δx = Δy. We choose Δx < dinj/2, where dinj is the distance between the centre of the injection borehole and the closest edge of the fault plane along strike. Here dinj ≈ 7 mm and we chose Δx = Δy = 2.5 mm. 3.2 Probabilistic inversion The probabilistic inversion is performed using the Metropolis Hastings algorithm [application of the Markov Chain Monte Carlo (MCMC) methods, Metropolis et al. 1953; Hastings 1970]. The concept of the method relies on conducting a large number of forward computations, with a new model for each computation. The acceptance probability for each new diffusivity model in estimated as: $$\begin{eqnarray*} A = \text{min} \left( \frac{{d}^{n+1}}{{d}^{n}}, 1 \right), \end{eqnarray*}$$(8) where d is the probability density function and the superscript n refers to number of the iteration. For the problem studied here: $$\begin{eqnarray*} {d} = \exp \left( - \frac{1}{2} \Sigma \left( \frac{{p}_\mathrm{obs}^N({t}) - {p}_\mathrm{obs}^E({t})}{\sigma ^\mathrm{obs}} \right)^2 \right), \end{eqnarray*}$$(9) where σobs is the uncertainty on the observational measurements in the laboratory. For this study, the uncertainty of the pressure sensors was σobs = 10−3 MPa. As the computation of the diffusion equation requires a rather small time step to assure numerical stability (Δt ≤ Δx2/2 max(D)), the diffusivity vector is parametrized on a very fine temporal grid. However, it is very complicated to apply the MCMC algorithm to such large number of parameters. For this reason, we apply it on a sparse diffusivity vector (formed by 6 model parameters here), that we later interpolate into the fine grid using a linear interpolation. We stop the inversion algorithm when the mean and the standard deviation of the accepted models do not evolve anymore. The acceptance rate of this inversion process is the ratio between the selected and the total models: nSel/nTot, and should be equal to 30 per cent. From the series of simulations, we can apply statistical analysis for all the accepted models, and display the distribution of each model parameter or the associated momentum and quantiles. As we expect an initial burn-in phase, in which the objective function rapidly decreases, until the model converges, all the statistics presented in the next section are estimated after this burn-in phase. 4 APPLICATION TO THE EXPERIMENTAL DATA We present here the results of the numerical inversion methods on the three experimental injection tests (Table 2). We start by exposing the application of the deterministic approach (adjoint state method) in order to estimate, in the least-squares sense, the best diffusivity model that can explain the experimental data, and then present the probabilistic approach (MCMC) so to estimate the uncertainty and the validity of the best model. We do not discuss the physical interpretation of the evolution of the hydraulic diffusivity here and we leave it to Section 5. 4.1 Estimating the best model: deterministic approach Figs 3, 4 and 5 illustrate the application of the deterministic approach (fit of the data and inverted diffusivity history) to the injection tests at 30, 60 and 95 MPa confining pressure, respectively. Additionally, in Appendix B are presented some additional results, in particular the objective functions and the reconstructed model of pore pressure. The number of non-linear iterations was set to 1400 at 30 MPa verifying that the relative reduction of the objective function was smaller than 10−3 (Fig. A2a). At 95 MPa, the number was set to 1500 iterations, in which the objective function decreased from 400 to 2×10−1 (Fig. A5a). As for the test at 60 MPa, results are presented at 200, 1000, 3000 and 5000 iterations. The initial diffusivity model was considered constant and time-independent, that is independent of the normal stress reduction or slip accumulation throughout the experiment, and was taken D0(t) = 4.6×10−6 m2 s–1 at 30 MPa, D0(t) = 1.6×10−6 m2 s–1 at 60 MPa, and D0(t) = 6.6×10−7 m2 s–1 at 95 MPa (Figs 3b, 4b and 5b). Several initial models were tested, without any important effect observed on the inversion results. Figure 3. Open in new tabDownload slide Results of the deterministic approach for the experiment at Pc = 30 MPa: (a) Numerical and experimental pressures at the injection and observation boreholes: black colour refers to the the experimental data, red colour to the numerical results, solid line represents the injection borehole and dashed ones represent the observation borehole. Areas 1, 2 and 3 are zoomed and represented in Fig. A3; (b) Time-series of the diffusivity model: the black line represents the best model (inversion result), and the black dashed line represents the initial model. Figure 3. Open in new tabDownload slide Results of the deterministic approach for the experiment at Pc = 30 MPa: (a) Numerical and experimental pressures at the injection and observation boreholes: black colour refers to the the experimental data, red colour to the numerical results, solid line represents the injection borehole and dashed ones represent the observation borehole. Areas 1, 2 and 3 are zoomed and represented in Fig. A3; (b) Time-series of the diffusivity model: the black line represents the best model (inversion result), and the black dashed line represents the initial model. Figure 4. Open in new tabDownload slide Same as for Fig. 3 for the injection test at 60 MPa confining pressure, except the coloured plots refer to the numerical results: with green, red, dark blue and light blue for 200, 1000, 3000 and 5000 inversion iterations, respectively. Figure 4. Open in new tabDownload slide Same as for Fig. 3 for the injection test at 60 MPa confining pressure, except the coloured plots refer to the numerical results: with green, red, dark blue and light blue for 200, 1000, 3000 and 5000 inversion iterations, respectively. Figure 5. Open in new tabDownload slide Same as for Fig. 3 for the injection test at 95 MPa confining pressure. Figure 5. Open in new tabDownload slide Same as for Fig. 3 for the injection test at 95 MPa confining pressure. The modelled pressure at the observation borehole p|$_\mathrm{obs}^N$| (red plot in Figs 3a, 4a and 5a) replicates relatively well the experimental measurement for the majority of the time domain, however some discrepancies are observed. We particularly comment the injection test at 30 MPa as the behaviour is quite similar for the different tests. Fig. A3 is a zoom version of Fig. 3 a over different time intervals. In time interval (1) (Fig. A3a): 180 < t < 280 s, at this time range, poro-elastic effects are suspected. The model underestimates the pressure, and we are not able to explain perfectly the experimental data using only pure diffusion. From Fig. 3(c), we observe in this time range an artificial bump in the best model of diffusivity, which can not be physically interpreted, as it may be an artificial way for the model to try and replicate the experimental data. For the time interval (2) (Fig. A3b), the measurement and the modelled pressure are quite similar. The time interval (3) (Fig. A3c) represents the largest times in the experiment t > 850 s, here also the model underestimates the measured pressure. It is difficult with our current model to fit the observations as the limited time remaining from the experiment may be not enough to properly model the diffusion process. From Fig. 3(c), we observe a plateau in the best model in this time range. We can not perfectly rely on the solution in this domain as the model was not able to fit perfectly the data here, however the plateau value can give us an estimation about the value of the hydraulic diffusivity in this time range. The experiment at 60 MPa is quite particular as it exhibits different step-like increases in pressure (Fig. 2b). We were not able to reproduce such a behaviour via modeling only a pure pressure diffusion processes. For this inversion, we show results after 200, 1000, 3000 and 5000 iterations. From Fig. 4(a), we can see that the numerical model fails to replicate the pressure steps, for instance at t ≈ 1250 and 1370 s. The diffusivity model after 200 iterations is rather smooth, whereas various peaks in the diffusivity model appear at the different locations of the pressure steps (Fig. 4b). The amplitudes of the peaks can be more important as the number of inversion iterations increases. We are not sure whether these peaks have a physical meaning and to what extent we can be confident as to their amplitude, the problem is underdetermined at these time intervals. We recall that no regularization of the diffusivity model has been applied. For these reasons, we choose to stop the inversion at 1000 iterations, even though the objective function is not at its minimum (Fig. A4a), so as to keep a relative smooth solution. For the injection test at 95 MPa, pump A was emptied at t ≈ 2500 s, causing the pressure at the injection borehole to suddenly drop to 0, then re-increase again when the pump was replaced (Fig. 5a). For this reason, we observe small artificial abrupt changes in the diffusivity around 2500 s (Fig. 5b). Here we chose to take into consideration the emptying of the pump so not to affect the diffusivity history. However, we could apply a mask over this time range to remove this part of the data, or apply regularization to the modeling parameters. For all the reasons mentioned before, we will only focus on the following time ranges when physically interpreting the evolution of the hydraulic diffusivity throughout the injection experiment: [≈ 300–800] s at 30 MPa, [700–1400] s at 60 MPa and [2000–4200] s at 95 MPa. In these time ranges, we can observe that the hydraulic diffusivity increases throughout the injection experiment, as the pressure diffuses along the fault plane. 4.2 Estimating the uncertainties: the MCMC approach The application of the MCMC method allows us to estimate to what extent the best model is valid and to explore the range of uncertainty for diffusivity history. We chose 6 model parameters to describe D(t), evenly spaced in time, and we started the parametrization in the time range where the diffusivity solution, issued from the deterministic approach, is not constant. For instance, the first model parameter is at t = 280 s for the test at 30 MPa (Fig. 6a). The initial model Dinit is taken constant, equal to the diffusivity value estimated by the deterministic approach before injection started. We then allow the algorithm to generate new models in the range [0.01–15]Dinit. We applied the MCMC method for 40 000 iterations, after which we verified that the average and the standard deviations, over the accepted models, of the 6 model parameters have reached stabilization. Fig. 6 illustrates a statistical representation the accepted models, in particular the probability density of the accepted models for the three injection tests at 30, 60 and 95 MPa. As mentioned in Section 3.2, the statistics are conducted for the accepted models after the burn-in phase, which is fixed at 100 accepted models. For the different injection tests, the deterministic best model is close the average |$\overline{D}$|(MCMC) for the majority of the time range: exceptions are observed at 60 MPa. We also note that the standard deviations are relatively small, and the quantiles are narrow, in the time intervals mentioned above. As expected, the range of accepted models is wider for the largest time of the different tests, as the model is not constrained is this part. We should also note that we do not expect to find the peaks and oscillations in the diffusivity model for the case at 60 MPa, due to the linear interpolation and as we did not impose inversion points at these locations. Figure 6. Open in new tabDownload slide Results of the MCMC Method for the experiment at Pc = (a) 30 MPa, (b) 60 MPa and (c) 95 MPa; Initial diffusivity model: red, average of the accepted models: grey (scattered stars represent the inverted points of MCMC, standard deviation: dashed grey lines, quantiles at 68 and 95 per cent: dashed green and pink lines, and the best model obtained with the deterministic method: light blue. The colour scale illustrates the probability density of the accepted models. Figure 6. Open in new tabDownload slide Results of the MCMC Method for the experiment at Pc = (a) 30 MPa, (b) 60 MPa and (c) 95 MPa; Initial diffusivity model: red, average of the accepted models: grey (scattered stars represent the inverted points of MCMC, standard deviation: dashed grey lines, quantiles at 68 and 95 per cent: dashed green and pink lines, and the best model obtained with the deterministic method: light blue. The colour scale illustrates the probability density of the accepted models. 4.3 Discussion We presented a deterministic inversion using the simple steepest descent algorithm to minimize the objective function, in which around 1000–1500 iterations were needed to converge, for which the computation time was quite acceptable. As it mainly depends on the resolution of the diffusion problem, it can be large if the diffusion problem is more complex. In this case, the application of different approaches, in which only a few dozens iterations would be needed can be quite effective, for instance: the conjugate gradient method (Shewchuk 1994), the quasi-Newton approaches (Kelley 1999), or the use of a pre-conditioner (Brossier 2011; Metivier et al. 2013). Here, no regularization was introduced to the deterministic inversion algorithm. Regularization techniques rely on additional information in order to stabilize the solution, for instance if the solution is supposed to be bound between some limits, smooth, have non-zero elements, etc. Different regularization techniques exist (Thikonov & Arsenin 1977; Engl et al. 2000; Kaltenbacher et al. 2008; Whitney 2009). More research is needed to determine the optimal regularization term as well as the values of the hyperparameters. The use of such methods might allow to avoid the numerical oscillations observed in the diffusivity history for the injection at 60 MPa confining pressure (Fig. 4b). Furthermore, we consider here a spatial homogeneous hydraulic diffusivity along the fault plane, as we only have one point of observation. As, the inversion method can be applied to multiple observation points, if the experimental protocol can be equipped by multiple pressure sensors on the fault plane, it would allow for the construction a 3-D (2-D space and time) diffusivity matrix. Here we explored the uncertainties with respect to the best model of the diffusivity using of a probabilistic inversion, in which we considered model parameters evenly spaced in time. As the variation of the hydraulic diffusivity coefficient is not the same through the injection experiment (Figs 3a, 5a and 4a), a variable spacing parametrization (in particular fine parametrization over the time intervals where the diffusivity varies a lot) may better reflect the uncertainties relative to the diffusivity variation. Another development is the update of the interpolation procedure between the model parameters from linear to a quadratic or Lagrangian interpolation. The method presented here can be applied to reservoir monitoring, and may be used to estimate hydraulic diffusivity variations. In this case, multiple measures of pressure are needed in addition to the one measured at the injection borehole, in order to better constrain the inverted diffusivity history. The inversion algorithm remains quite the same. The additional feature to the presented model is the need to implement a more complex fluid flow model in order to consider a bulk permeability, the presence of different fluids and their possible interaction, and eventually the pressure diffusion process into a fault network. In the next section, we investigate the hydraulic diffusivity variation in the time intervals mentioned above. We use as the diffusivity time history the best model issued from the deterministic approach, and consider the standard deviations calculated using the MCMC method as an estimation of the uncertainties in the model. 5 DIFFUSIVITY, DISPLACEMENT AND EFFECTIVE STRESS The estimated hydraulic diffusivity varies in the ranges [4.6×10−6 – 3.7×10−5], [1.6×10−6 – 1.9×10−5] and [6.6×10−7 – 5.3×10−6] m2 s–1, at 30, 60 and 95 MPa confining pressure, respectively. A dependence on confining pressure emerges, as we observe systematically lower values at greater confining pressures. This was previously observed in different studies (e.g., Zoback & Byerlee 1975; Wibberley & Shimamoto 2003). Through the injection experiment, the hydraulic diffusivity increases by one order of magnitude approximately, as shear displacement accumulates on the fault and effective stress decreases following pressure diffusion. Following Jaeger et al. (2007), such diffusivity changes could be interpreted as permeability enhancement. Various observations associate permeability enhancement with either effective stress reduction (e.g. Zoback & Byerlee 1975; McKee et al. 1988; Fisher & Zwart 1996; Ghabezloo et al. 2009), or with slip accumulation (e.g. Esaki et al. 1999; Baghbanan & Jing 2008; Guglielmi et al. 2015b; Bhattacharya & Viesca 2019). Here we raise the question on how these two parameters could influence diffusivity enhancement in the same context. Fig. 7 illustrates the relation between the hydraulic diffusivity variation and the mean effective stress which is estimated as σeff = σn − pm, where σn is the average normal stress along the fault plane estimated using eq. (1) and pm is the mean pressure along the fault plane, estimated from the reconstructed pressure profile, using the best model issued from the deterministic method (Figs A2b, A4b and A5b). The results from the three injection tests present the same dependence and tendency: a decrease of the hydraulic diffusivity with the effective normal stress. At 30 and 95 MPa confining pressure, the results collapse relatively on the same curve, while the ones at 60 MPa present a small difference. This dependence of the hydraulic properties of the rock to the effective stress may be affected by the opening/closure of the joints along the fracture plane. According to Bandis et al. (1983), the effective normal stress and the joint closure are related as follows: $$\begin{eqnarray*} \ln \frac{\sigma }{\sigma _0} = \mathrm{J} \delta , \end{eqnarray*}$$(10) where σ0 represents the initial normal stress, J is a constant with dimensions of inverse length and δ is the joint closure. If we (1) assume that δ = h – h0, where h is the hydraulic aperture and h0 is the initial aperture, (2) consider the cubic law where K = h2/12 (K being the permeability, Snow 1965; Witherspoon et al. 1980; Bradley 1987) and (3) consider the following relation between the permeability and the hydraulic diffusivity D = k/(ϕμc), where ϕ is the porosity, μ is the viscosity and c is the total compressibility (Jaeger et al. 2007), we get the following relation between the hydraulic diffusivity and the effective stress $$\begin{eqnarray*} \frac{{D}}{{D}_0} = \left( 1 - \frac{1}{{J} {h}_0} \ln \frac{\sigma }{\sigma _0}\right)^2, \end{eqnarray*}$$(11) where D0 represents the initial diffusivity value, before the injection test started. Fig. 8 represents the relation between the square-root of D/D0 and the logarithm of the ratio σ/σ0. Hence, the dependence between the hydraulic diffusivity and the effective normal stress appears to follow the previous relation (eq. 11). At 30 and 95 MPa of confining pressure, the results reflect more or less the same behaviour (1/Jh0 ≈ 1.5). A clear difference is however observed at 60 MPa (1/Jh0 ≈ 2.1). This observation suggests that not only the effective stress can affect the hydraulic diffusivity, an additional factor might have an impact as well. Figure 7. Open in new tabDownload slide Hydraulic diffusivity and mean effective stress in the time ranges: [300–700] s at 30 MPa, [700–1400] s at 60 MPa, and [2000–4200] s at 95 MPa. The effective stress is computed as the difference between the normal stress and the mean pore pressure along the fault which is numerically estimated using the inverted diffusivity history. Different colours refer to the different confining pressure values. Diffusivity values are issued from the deterministic method, and the error-bars represented here are the standard deviations estimated using the MCMC method. Figure 7. Open in new tabDownload slide Hydraulic diffusivity and mean effective stress in the time ranges: [300–700] s at 30 MPa, [700–1400] s at 60 MPa, and [2000–4200] s at 95 MPa. The effective stress is computed as the difference between the normal stress and the mean pore pressure along the fault which is numerically estimated using the inverted diffusivity history. Different colours refer to the different confining pressure values. Diffusivity values are issued from the deterministic method, and the error-bars represented here are the standard deviations estimated using the MCMC method. Figure 8. Open in new tabDownload slide Squares of the ratio of the hydraulic diffusivity to the initial diffusivity and the logarithm of the ratio of the mean effective stress to the initial effective stress. Legend same as Fig. 7. Figure 8. Open in new tabDownload slide Squares of the ratio of the hydraulic diffusivity to the initial diffusivity and the logarithm of the ratio of the mean effective stress to the initial effective stress. Legend same as Fig. 7. According to Chen et al. (2000b), the hydraulic aperture depends on the shear displacement and the confining pressure. Fig. 9 illustrates the hydraulic diffusivity variation with average slip along the fault (only the injection phase is considered). The diffusivity increase is not directly dependent on the amount of slip: more slip for 95 MPa does not induce a larger diffusivity increase. For the injection tests at 30 and 95 MPa, the hydraulic diffusivity increases systematically with shear slip. However, this is not the case for the test at 60 MPa, as this experiment exhibits various stick-slip events as well as time periods of negligible displacement (Fig. 2b). During the latter, the diffusivity continues to increase with the reduction of effective stress independently of the shear slip, while during the stick slip events, an acceleration in the diffusivity increase is observed, that is later recompensated by a smaller decrease. We are however not sure to what extent we can rely on the quantitative abrupt increase in the hydraulic diffusivity, as discussed in the previous section. It appears that the nature of the slip events may affect the hydraulic diffusivity evolution, however more extensive work and future studies are needed to investigate this feature. In the context of a decametric scale injection test in the underground LSBB laboratory in the South France, Guglielmi et al. (2015a) also reported a difference in the permeability enhancement between the seismic and aseismic slip phases. Figure 9. Open in new tabDownload slide Hydraulic diffusivity variation and average slip. Legend same as Fig. 7. Figure 9. Open in new tabDownload slide Hydraulic diffusivity variation and average slip. Legend same as Fig. 7. Our investigation of the relation between the shear slip and diffusivity enhancement remains however limited, as we can only estimate an average shear slip along the fault plane and only invert for an effective diffusivity. Following stick slip events, localized shear slip along the fault plane might be expected, moreover, permeability anisotropy has been in fact observed to be induced by shear displacement (Zhang & Tullis 1998; Auradou et al. 2005). For this reason, we strongly recommend for future laboratory experiments to equip the rock sample with multiple pressure and displacement sensors. Such data would allow for the investigation of the hydraulic diffusivity enhancement accompanying stick slip events. The injection tests presented in this study are conducted for a saw-cut fault in a Andesite rock sample. Thus the results presented here correspond to a smooth surface. According to Ye & Ghassemi (2018), shear slip and permeability enhancement are largely dependent on the roughness of the fault. It would be interesting for future studies to investigate hydraulic diffusivity variations on experimental faults with different surface roughness. In a similar context of laboratory injection experiments, Passelegue et al. (2018) showed that the injection pressure rate controls the fluid pressures perturbation and can affect the fault reactivation. From a numerical perspective, Almakari et al. (2019) reported the important effect of the injection pressure rate on the seismicity rate and the magnitude content of the seismic events. In this study, fluid was injected into the fault under a constant injection pressure rate. Various injection pressure rates should be tested in future studies to investigate if the permeability enhancement may be affected as well. Finally, in addition to the characterization of the evolution of the hydraulic diffusivity through a single injection experiment and the investigation of how it varies with slip accumulation and effective stress reduction, the application of such numerical approaches could be very advantageous, as it potentially allows for reconstruction of the spatio-temporal pore pressure changes during the injection test, that allows for the track the diffusive front. With the use of appropriate acoustic sensors, this can give more insights into the relation between the diffusive front and the aseismic/seismic rupture front. 6 CONCLUSION We used in the context of this study laboratory experimental constant-rate injection tests on a saw-cut fault in an Andesite rock sample, under triaxial loading, at different values of confining pressure. The experimental sample was equipped by two pressure sensors that allow for the measurements of the pressure history at two locations along the fault plane. The fluid injection reactivated the experimental fault, inducing different slip events, that were accompanied by pressure drops in the injection borehole and pressure steps in the observation borehole, which could indicate an enhancement of the hydraulic diffusivity along the fault. We then developed and performed deterministic and probabilistic inversions to the experimental data, in particular the pore pressure history. This allows us to characterize the time history of an effective hydraulic diffusivity D(t) throughout the injection test. Deterministic inversion, using the adjoint state method, was performed to estimate the best diffusivity model that could explain the experimental data. Then, probabilistic inversion, using an MCMC algorithm, was applied in order to estimate the associated uncertainties. The numerical method was able to reproduce the experimental data for a wide time range of the different experiments. In the injection reactivation tests presented here, an enhancement of the hydraulic diffusivity by one order of magnitude was observed. Hydraulic diffusivity variations were found to strongly depend on the reduction of the normal effective stress acting along the fault plane, while a second-order effect of the shear slip was observed. More extensive work in future studies is needed to fully explore the effect of the nature of shear slip events on the hydraulic diffusivity evolution. Finally, the use of such method can be of broad interest at the reservoir scale in order to monitor hydraulic diffusivity variations. As well, with the use of ultrasonic-transducers, such method would allow to further investigate the relationship between fluid and fault reactivation. REFERENCES Almakari M. , Dublanchet P., Chauris H., Pellet F., 2019 . Effect of the injection-scenatio on the rate and magnitude content of injection-induced seismicity: case of a heterogeneous fault , J. geophys. Res. , 124 ( 8 ), 8426 – 8448 . 10.1029/2019JB017898 Google Scholar Crossref Search ADS WorldCat Auradou H. , Drazer G., Hulin J., Koplik J., 2005 . Permeability anisotropy induced by the shear displacement of rough fracture walls , Water Resour. Res. , 41 , W09423 , doi.org/10.1029/2005WR003938 . 10.1029/2005WR003938 Google Scholar Crossref Search ADS WorldCat Baghbanan A. , Jing L., 2008 . Stress effects on permeability in a fractured rock mass with correlated fracture length and aperture , Int. J. Rock Mech. Min. Sci. , 45 , 1320 – 1334 ., doi.org/10.1016/j.ijrmms.2008.01.015 . 10.1016/j.ijrmms.2008.01.015 Google Scholar Crossref Search ADS WorldCat Bandis S.C. , Lumsden A., Barton N.R., 1983 . Fundamentals of rock joint deformation , Int. J. Rock Mech. Min. Sci., Geomech. , 20 ( 6 ), 249 – 268 . Google Scholar Crossref Search ADS WorldCat Bhattacharya P. , Viesca R., 2019 . Fluid induced aseismic fault slip outpaces pore fluid migration , Science , 364 , 464 – 468 . Google Scholar Crossref Search ADS PubMed WorldCat Bradley H. , 1987 . Petroleum Engineering Handbook , Society of Petroleum Engineers . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Brossier R. , 2011 . Two-dimensional frequency-domain visco-elastic full waveform inversion: parallel algorithms, optimization and performance , Comput. Geosci. , 37 , 444 – 455 . Google Scholar Crossref Search ADS WorldCat Chavent G. , 1991 . On the theory and practice of non-linear least-squares , Adv. Water Resour. , 14 ( 2 ), 55 . Google Scholar Crossref Search ADS WorldCat Chen Z. , Narayan S., Yang Z., Rahman S., 2000a . An experimental investigation of hydraulic behaviour of fractures and joints in granitic rock , Int. J. Rock Mech. Min. Sci. , 37 , 1061 – 1071 . Google Scholar Crossref Search ADS WorldCat Chen Z. , Narayan S.P., Yang Z., Rahman S., 2000b . An experimental investigation of hydraulic behaviour of fractures and joints in granitic rock , Int. J. Rock Mech. Min. Sci. , 37 , 1061 – 1071 . Google Scholar Crossref Search ADS WorldCat Darcy H. , 1857 . Recherches experimentales relatives au mouveament de l’eau dans le tuyaux , Mallet - Bachelier, Quai des Augustins, 55, ecole imperiale polytechnique edn . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC De Barros L. , Guglielmi Y., Rivet D., Cappa F., Duboeuf L., 2018 . Seismicity and fault aseismic deformation caused by fluid injectionin decametric in-situ experiments , Compte Rend. Geosci. , 350 ( 8) , 464 – 475 . Google Scholar Crossref Search ADS WorldCat Dublanchet P. , 2019 . Fluid driven shear cracks on a strengthening rate-and-state frictional fault , J. Mech. Phys. Solids , 132 , 103672 . Google Scholar Crossref Search ADS WorldCat Duboeuf L. , De Barros L., Cappa F., Guglielmi Y., Deschamps A., Seguy S., 2017 . Aseismic motions drive a sparse seismicity during fluid injections into a fractured zone in a carbonate reservoir , J. geophy. Res. , 122 . Google Scholar OpenURL Placeholder Text WorldCat Engl H.W. , Hanke M., Neubauer A., 2000 . Regularization of Inverse Problems , Vol. 35 , Kluwer Academic Publishers . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Esaki T. , Du S., Mitani Y., Ikusuda K., Jing L., 1999 . Development of a shear-flow test apparatus and determination of coupled properties for a single rock joint , Int. J. Rock Mech. Min. Sci. , 36 ( 5 ), 641 – 650 . Google Scholar Crossref Search ADS WorldCat Fisher A. , Zwart G., 1996 . Relation between permeability and effective stress along a plate-boundary fault, Barbados accretionary complex , Geology , 24 ( 4 ), 307 – 310 . Google Scholar Crossref Search ADS WorldCat Galis M. , Ampuero J.-P., Mai P., Cappa F., 2017 . Induced seismicity provides insight into why earthquake ruptures stop , Sci. Adv. , 3 ( 12) , eaap7528, doi:10.1126/sciadv.aap7528 . Google Scholar OpenURL Placeholder Text WorldCat Garagash D. , Germanovish L., 2012 . Nucleation and arrest of dynamic slip on a pressurized fault , J. geophys. Res. , 117 ( B10310 ), doi:10.1029/2012JB009209. Google Scholar OpenURL Placeholder Text WorldCat Ghabezloo S. , Sulem J., Guedon S., Martineau F., 2009 . Effective stress law for the permeability of a limestone , Int. J. Rock Mech. Min. Sci. , 46 , 297 – 306 . Google Scholar Crossref Search ADS WorldCat Guglielmi Y. , Cappa F., Avouac J.-P., Henry P., Elsworth D., 2015a . Seismicity triggered by fluid injection induced aseismic slip , Science , 348 , 1224 – 1226 . Google Scholar Crossref Search ADS WorldCat Guglielmi Y. , Elsworth D., Cappa F., Henry P., Gout C., Dick P., Durand J., 2015b . In situ observations on the coupling between hydraulic diffusivity and displacements during fault reactivation in shales , J. geophys. Res. , 120 , 7729 – 7748 . Google Scholar Crossref Search ADS WorldCat Gutierrez M. , øino L., Nygard R., 2000 . Stress-dependent permeability of a de-mineralised fracture in shale , Mar. Petrol. Geol. , 17 , 895 – 907 . Google Scholar Crossref Search ADS WorldCat Hastings W. , 1970 . Monte Carlo sampling methods using Markov chains and their applications , Biometrika , 57 ( 1 ), 97 . Google Scholar Crossref Search ADS WorldCat Healy J. , Rubey W., Griggs D., Raleigh C., 1968 . Disposal of waste fluids by injection into a deep well has triggered earthquakes near Denver, Colorado , Science , 161 ( 3848 ), 1301 – 1310 . Google Scholar Crossref Search ADS PubMed WorldCat Hermann R. , Park S.-K., Wang C.-Y., 1981 . The Denver earthquakes of 1967-1968 , Bull. seism. Soc. Am. , 71 ( 3 ), 731 – 745 . Google Scholar OpenURL Placeholder Text WorldCat Hsieh P. , Bredehoeft J., 1981 . A reservoir analysis of the Denver earthquakes: a case of induced seismicity , J. geophys. Res. , 86 ( No. B2 ), 903 – 920 . Google Scholar Crossref Search ADS WorldCat Im K. , Elsworth D., Fang Y., 2018 . The influence of preslip sealing on the permeability evolution of fractures and faults , Geophys. Res. Lett. , 45 , 166 – 175 . Google Scholar Crossref Search ADS WorldCat Jaeger J. , Cook N., Zimmerman R., 2007 . Fundamentals of Rock Mechanics , Blackwell Publishing . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Jarny Y. , Ozisik M., Bardon J., 1991 . A general optimization method using adjoint equation for solving multidimensional inverse heat conduction , Int. J. Heat Mass Transfer , 134 ( 11 ), 2911 – 2919 . Google Scholar Crossref Search ADS WorldCat Kaltenbacher B. , Neubauer A., Scherzer O., 2008 . Iterative Resularization Methods for Nonlinear Problems , De Gruyter . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Kelley C. , 1999 . Iterative Methods for Optimization , SIAM . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Kim W.-Y. , 2013 . Induced seismicity associated with fluid injection into a deep well in Youngstown, Ohio , J. geophys. Res. , 118 , 3506 – 3518 . Google Scholar Crossref Search ADS WorldCat Li Z. , Fortin F., Nicolas A., Deldicque D., Guguen Y., 2019 . Physical and mechanical properties of thermally cracked andesite under pressure , Rock Mech. Rock Eng. , 52 , 3509 – 3529 . Google Scholar Crossref Search ADS WorldCat McKee C. , Bumb A., Koelng R., 1988 . Stress-dependent permeability and porosity of coal and other geologic formations , SPE Format. Eval. , 3 ( 01) , doi:10.2118/12858-PA . Google Scholar OpenURL Placeholder Text WorldCat Metivier L. , Brossier R., Virieux J., Operto S., 2013 . Full waveform inversion and the truncated Newton method , SIAM J. Scient. Comput. , 35 ( 2 ). Google Scholar OpenURL Placeholder Text WorldCat Metropolis N. , Rosenbluth A., Rosenbluth M., Telles A., 1953 . Equation of state calculations by fast computing machines , J. Chem. Phys. , 21 , 1087 . Google Scholar Crossref Search ADS WorldCat Passelegue F. , Brantut N., Mitchell T., 2018 . Fault reactivation by fluid injection: controls from stress state and injection rate , Geophys. Res. Lett. , 45 , 12837 – 12846 . Google Scholar OpenURL Placeholder Text WorldCat Plessix R.-E. , 2006 . A review of the adjoint-state method for computing the gradient of a functional with geophysical applications , Geophys. J. Int. , 167 , 495 – 503 . Google Scholar Crossref Search ADS WorldCat Robert C. , Casilla G., 2004 . Monte Carlo Statistical Methods , Springer . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Rutter E. , Mecklenburgh J., 2018 . Influence of normal and shear stress on the hydraulic transmissivity of thin cracks in a tight quartz sandstone, a granite, and a shale , J. geophys. Res. , 123 , 1262 – 1285 . Google Scholar Crossref Search ADS WorldCat Shapiro S. , Huenges E., Borm G., 1997 . Estimating the crust permeability from fluid-injection-induced seismic emission at the KTB site , Geophys. J. Int. , 131 ( 2 ), F15 – F18 . Google Scholar Crossref Search ADS WorldCat Shapiro S. , Rothert E., Rath V., Rindschwentner J., 2002 . Characterization of fluid transport properties of reservoirs using induced microseismicity , Geophysics , 67 ( 1 ), 212 – 220 . Google Scholar Crossref Search ADS WorldCat Shewchuk J. , 1994 . An Introduction to the Conjugate Gradient Method Without the Agonizing Pain , Carnegie Mellon University . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Snow D. , 1965 , A parallel plate model of fractures permeable media , PhD thesis , University of California , Berkeley . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Thikonov A. , Arsenin V., 1977 . Solution of Ill-Posed Problems , Winston . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Whitney M.L. , 2009 , Theoretical and numerical study of Tikonov’s regularization and Morozov’s discrepancy principle . PhD thesis , Georgia State University, Department of Mathematics and Statistics . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wibberley C. , Shimamoto T., 2003 . Internal structure and permeability of major strike-slip fault zones: the median tectonic line in mie prefecture, Southwest Japan , J. Struct. Geol. , 25 , 59 – 78 . Google Scholar Crossref Search ADS WorldCat Witherspoon P.A. , Wang J.S.Y., Iwai K., Gale J.E., 1980 . Validity of cubic law for fluid in a deformable rock fracture , Water Resour. Res. , 16 ( 6 ), 1016 – 1024 . Google Scholar Crossref Search ADS WorldCat Wu W. , Reece J., Gensterblum Y., Zoback M., 2017 . Permeability evolution of slowly slipping faults in shale reservoirs , Geophys. Res. Lett. , 44 , 11368 – 11375 . Google Scholar Crossref Search ADS WorldCat Ye Z. , Ghassemi A., 2018 . Injection-Induced shear slip and permeability enhancement in granite fractures , J. geophy. Res. , 123 , 9009 – 9032 . Google Scholar Crossref Search ADS WorldCat Zhang S. , Tullis T., 1998 . The effect of fault slip on permeability and permeability anisotropy in quartz gouge , Tectonophysics , 295 , 41 – 52 . Google Scholar Crossref Search ADS WorldCat Zoback M.D. , Byerlee J.D., 1975 . Permeability and effective stress , Geologic Notes , 59 ( 1 ), 154 – 158 . Google Scholar OpenURL Placeholder Text WorldCat APPENDIX A: DEVELOPMENT OF THE ADJOINT STATE METHOD We present here the development of the adjoint state method applied to the inversion of the hydraulic diffusivity history in the context of a 2-D pressure diffusion inside an elliptical fault plane. We use this method to estimate the gradient of the objective function. In the context of this study, pressure p|$_\mathrm{inj}^E$|(t) is injected in the injection borehole (coordinates xinj,yinj), and measured p|$_\mathrm{obs}^E$|(t) in the observation borehole (coordinates xobs,yobs). We consider an effective diffusivity D(t) along the fault plane, which is estimated in a least-squares sense (between the experimental value of the pore pressure at the observation borehole and the numerical one, Chavent 1991; Jarny et al. 1991). It means we assume an additive noise with a Gaussian distribution. The objective function reads: $$\begin{eqnarray*} J[D({t})] = \frac{1}{2} \int \left( {p}_\mathrm{obs}^N[D]({t}) - {p}_\mathrm{obs}^E({t}) \right)^2 \mathrm{d}t. \end{eqnarray*}$$(A1) The adjoint state method introduces a new variable λ in the estimation of the objective function, defined as: $$\begin{eqnarray*} J[D,\lambda ,{p}] &=& \frac{1}{2} \iiint \left( \delta \left( {x} - {x}_\mathrm{obs} \right) \delta \left( {y} - {y}_\mathrm{obs} \right) \left( {p}_\mathrm{obs}^N - {p}_\mathrm{obs}^E \right)^2 \right) {{\rm d}x~{\rm d}y~{\rm d}t} \\ && - \, \iiint \lambda ({x},{y},{t}) . \left[ \frac{\partial {p}}{\partial {t}} - D \Delta {p} \right] {{\rm d}x~{\rm d}y~{\rm d}t}. \end{eqnarray*}$$(A2) For the sake of simplicity, we did not include the initial and boundary conditions in the previous equation. The total gradient of the objective function is: $$\begin{eqnarray*} \frac{{{\rm d}J}}{{{\rm d}D}} = \frac{\partial J}{\partial D} + \frac{\partial \lambda }{\partial D} \frac{\partial J}{\partial \lambda } + \frac{\partial \mathrm{p}}{\partial D} \frac{\partial J}{\partial {p}}. \end{eqnarray*}$$(A3) where D = D(t), p = p(x,y,t) and λ = λ(x,y,t). We choose p and λ so that ∂J/∂p = 0 and ∂J/∂λ = 0, to avoid estimating the Fréchet derivatives ∂p/∂D and ∂λ/∂D. Using eq. (A2): $$\begin{eqnarray*} \frac{\partial J}{\partial D({t_0})} = \iint \lambda ({x},{y},{t}_0) \Delta {p}({x},{y},{t}_0) {{\rm d}x~{\rm d}y}. \end{eqnarray*}$$(A4) Thus the gradient of the objective function becomes: $$\begin{eqnarray*} \frac{{{\rm d}J}}{{{\rm d}D}} = \iint \lambda ~\Delta {p}~{{\rm d}x~{\rm d}y}. \end{eqnarray*}$$(A5) In order to estimate the gradient, we need the pressure p and the adjoint state variable λ. We thus estimate the partial derivatives of the objective function with respect to p and λ using eq. (A2): $$\begin{eqnarray*} \frac{\partial J}{\partial \lambda ({x_0,y_0,t_0})} = - \left( \frac{\partial {p}(\mathrm{x_0,y_0,t_0})}{\partial {t}} - D({t}_0) \Delta {p}({x_0,y_0,t_0}) \right). \end{eqnarray*}$$(A6) As we impose ∂J/∂λ = 0, thus eq. (A6) becomes: $$\begin{eqnarray*} \frac{\partial {p}}{\partial {t}} - D({t}) \Delta {p} = 0, \end{eqnarray*}$$(A7) which is the pressure diffusion equation, the state equation (by definition). For the derivative with respect to p, we first need to integrate by part. The partial derivative writes: $$\begin{eqnarray*} && \frac{\partial J}{\partial {p}({x,y,t})} = \\ && \frac{\partial }{\partial \text{p}({x,y,t})} \frac{1}{2} \iiint \Big ( \delta ({x - x_{\rm obs}}) \delta ({y - y_{\rm obs}}) \left( {p_{\rm obs}^N} - {p_{\rm obs}^E} \right)\Big )^2 {{\rm d}x~{\rm d}y~{\rm d}t} \\ && - \, \frac{\partial }{\partial {p}({x_0,y_0,t_0})} \iiint \lambda ({x,y,t}) \left[ \frac{\partial {p}}{\partial {t}} \right] {{\rm d}x~{\rm d}y~{\rm d}t} \\ && + \, \frac{\partial }{\partial {p}({x_0,y_0,t_0})} \iiint \lambda ({x,y,t}) D \left( \frac{\partial ^2 {p}}{\partial {x}^2} + \frac{\partial ^2 {p}}{\partial {y}^2} \right) {{\rm d}x~{\rm d}y~{\rm d}t} . \end{eqnarray*}$$(A8) We will now resolve independently each term of the right part of eq. (A8). For the sake of simplicity, we will name the 3 parts A, B and C, that is ∂J/∂p(x, y, t) = A – B + C. Part 1 $$\begin{eqnarray*} A = \delta ({x - x_{\rm obs}}) \delta ({y - y_{\rm obs}}) \left( {p_{\rm obs}^N} - {p_{\rm obs}^E} \right). \end{eqnarray*}$$(A9) Part 2: To solve this part we conduct an integration by part with respect to time t. $$\begin{eqnarray*} B &=& \frac{\partial }{\partial {p}({x_0,y_0,t_0})} \iint \Big [ \lambda ({x,y,t}) {p}({x,y,t}) \Big ]_0^{T_{\mathrm{max}}} {{\rm d}x~{\rm d}y} \\ && - \, \frac{\partial }{\partial {p}({x_0,y_0,t_0})} \iiint \frac{\partial \lambda }{\partial {t}} {p}({x,y,t}) {{\rm d}x~{\rm d}y~{\rm d}t} , \end{eqnarray*}$$(A10) where Tmax is the time at the end of the experiment. The pressure p is 0 at the start of the injection (p(x,y,t0) = 0), and we impose that λ at t = Tmax is zero (λ(x,y,Tmax) = 0), thus the first term of the right-hand part of the previous equation is 0. Thus, $$\begin{eqnarray*} B = - \frac{\partial }{\partial {p}({x_0,y_0,t_0})} \iiint \frac{\partial \lambda }{\partial {t}} {p}({x,y,t}) {{\rm d}x~{\rm d}y~{\rm d}t} = \frac{\partial \lambda ({x_0,y_0,t_0})}{\partial {t}} \end{eqnarray*}$$(A11) Part 3: To solve this part we use Green’s identities. We know that: $$\begin{eqnarray*} \frac{\partial }{\partial {x}} \left( \lambda \frac{\partial {p}}{\partial {x}} \right) = \frac{\partial \lambda }{\partial {x}} \frac{\partial {p}}{\partial {x}} + \lambda \frac{\partial ^2 {p}}{\partial {x}^2}, \end{eqnarray*}$$(A12) $$\begin{eqnarray*} \frac{\partial }{\partial {y}} \left( \lambda \frac{\partial {p}}{\partial {y}} \right) = \frac{\partial \lambda }{\partial {y}} \frac{\partial {p}}{\partial {y}} + \lambda \frac{\partial ^2 {p}}{\partial {y}^2}. \end{eqnarray*}$$(A13) By integrating eqs (A12) and (A13) over the domain Ω, we get: $$\begin{eqnarray*} \int _\Omega \overline{\div } \left( \lambda \overline{\nabla } {p} \right) {{\rm d}S} = \int _\Omega \overline{\nabla } \lambda . \overline{\nabla } {p} {{\rm d}S} + \int _\Omega \lambda \Delta {p} {{\rm d}S}, \end{eqnarray*}$$(A14) where dS = dx dy. This gives: $$\begin{eqnarray*} \int _{{{\rm d}} \Omega } \lambda ~\overline{\nabla } {p}~. \overline{{n}}~{{\rm d}u} = \int _\Omega \overline{\nabla } \lambda ~. \overline{\nabla } {p}~{{\rm d}S} + \int _\Omega \lambda ~\Delta {p}~{{\rm d}S}, \end{eqnarray*}$$(A15) where |$\overline{{n}}$| represents the normal vector with respect to the contour dΩ, and du is the distance along the contour dΩ. As no liquid is allowed to flow outside of the fault plane, the flow at the boundaries of the fault is 0 and Neuman boundary conditions are imposed: |$\overline{\nabla } {p} . \overline{{n}}$| = 0 along dΩ, giving: $$\begin{eqnarray*} \int _\Omega \overline{\nabla } \lambda ~. \overline{\nabla } {p}~{{\rm d}S} + \int _\Omega \lambda ~\Delta {p}~{{\rm d}S} = 0. \end{eqnarray*}$$(A16) Using Green’s identities to the adjoint state variable λ and applying the same reasoning, we get: $$\begin{eqnarray*} \int _{{\rm d} \Omega } {p}~\overline{\nabla } \lambda ~. \overline{{n}}~{{\rm d}u} = \int _\Omega \overline{\nabla } {p}~. \overline{\nabla } \lambda ~{{\rm d}S} + \int _\Omega {p}~\Delta \lambda ~{{\rm d}S}, \end{eqnarray*}$$(A17) We impose no flow for the adjoint state variable λ at the boundaries of the fault plane, thus |$\overline{\nabla } \lambda . \overline{{n}}$| = 0 along dΩ, giving: $$\begin{eqnarray*} \int _\Omega \overline{\nabla } {p}~. \overline{\nabla } \lambda ~{{\rm d}S} + \int _\Omega {p}~\Delta \lambda ~{{\rm d}S} = 0. \end{eqnarray*}$$(A18) By substracting eq. (A16) from eq. (A18), we get: $$\begin{eqnarray*} \int _\Omega \lambda ~\Delta {p}~{{\rm d}S} = \int _\Omega {p}~\Delta \lambda ~{{\rm d}S}. \end{eqnarray*}$$(A19) If we substitute eq. (A19) into term C, we get: $$\begin{eqnarray*} C &=& \frac{\partial }{\partial {p}({x_0,y_0,t_0})} \left( \iiint \left( \frac{\partial ^2 \lambda }{\partial {x}^2} + \frac{\partial ^2 \lambda }{\partial {y}^2} \right) D({t}) {p}({x,y,t}) {{\rm d}x~{\rm d}y~{\rm d}t} \right) \\ &=& \left( \frac{\partial ^2 \lambda }{\partial {x}^2} + \frac{\partial ^2 \lambda }{\partial {y}^2} \right) D({t}) \\ &=& \Delta \lambda ~D({t}). \end{eqnarray*}$$(A20) Therefore, assuming ∂J/∂p = 0, eq. (A8) becomes: $$\begin{eqnarray*} \frac{\partial \lambda }{\partial {t}} + D({t})~\Delta \lambda = \delta ({x - x_{\rm obs}}) \delta ({y - y_{\rm obs}}) \left( {p_{\rm obs}^N} - {p_{\rm obs}^E} \right) . \end{eqnarray*}$$(A21) Eqs (A7) and (A21) form a set of differential equations for the pressure p and the adjoint state variable λ. Once resolved, their solutions can be implemented in eq. (A5) to estimate the gradient of the objective function. A1 Validation of the gradient estimation The final step is to validate this method for the gradient estimation. We first need to verify the key parameters, that are the initial and boundary conditions as well as the source terms for the pressure p and the adjoint state variable λ. Then we estimate the gradient using the finite difference method, as follows: $$\begin{eqnarray*} \frac{\partial J}{\partial D{(t)}} = \frac{ J\left( D({t}) + \delta D({t}) \right) - J\left( D({t}) - \delta D({t}) \right) }{2 \delta D({t})}. \end{eqnarray*}$$(A22) To do so, we need to impose a local perturbation δD(t) at each time step t0, and then solve eq. (A7) once with the diffusivity model D(t)+δ D(t) and then D(t)–δ D(t). We estimate for each the residuals and the objective function using eq. (A1) and then implement the values in eq. (A22). This method is very long as it requires the resolution of the forward problem twice for each time step, for each inversion iteration. This is why we only use it for one inversion iteration, and for a few time steps in order to validate the gradient estimation by the adjoint state method. Fig. A1 illustrates a validation test that we conducted for a synthetic test: the gradients estimated using the adjoint state formulation and the finite difference method for a synthetic test are similar. Figure A1. Open in new tabDownload slide Validation of the gradient estimation for a synthetic test. In red: the normalized gradient estimated using the adjoint state formulation; Black circles: the normalized gradient estimated using the finite difference method. These values are for one inversion iteration. Figure A1. Open in new tabDownload slide Validation of the gradient estimation for a synthetic test. In red: the normalized gradient estimated using the adjoint state formulation; Black circles: the normalized gradient estimated using the finite difference method. These values are for one inversion iteration. Figure A2. Open in new tabDownload slide Detailed results of the deterministic approach for the experiment at Pc = 30 MPa: (a) Objective function on a semi-logarithmic scale. Stopping criteria: 1400 iterations and the verification that the reduction of the objective function (J(n) – J(n+1))/J(n) becomes smaller than 10−3; (b) Pore pressure profiles along strike of the fault in the reconstructed model: the colour scale refers to the time, a profile is plotted each ∼42 s. The vertical dashed lines represent the locations of the injection and observation boreholes, and the grey areas are the boundaries of the fault plane where no fluid flow is imposed. Figure A2. Open in new tabDownload slide Detailed results of the deterministic approach for the experiment at Pc = 30 MPa: (a) Objective function on a semi-logarithmic scale. Stopping criteria: 1400 iterations and the verification that the reduction of the objective function (J(n) – J(n+1))/J(n) becomes smaller than 10−3; (b) Pore pressure profiles along strike of the fault in the reconstructed model: the colour scale refers to the time, a profile is plotted each ∼42 s. The vertical dashed lines represent the locations of the injection and observation boreholes, and the grey areas are the boundaries of the fault plane where no fluid flow is imposed. Figure A3. Open in new tabDownload slide Zooms of Fig. 3(a) for only the pressures at the observation borehole, in the following time ranges: (a) 100 < t < 400 s; (b) 400 < t < 700 s; (c) 800 < t < 1100 s. Figure A3. Open in new tabDownload slide Zooms of Fig. 3(a) for only the pressures at the observation borehole, in the following time ranges: (a) 100 < t < 400 s; (b) 400 < t < 700 s; (c) 800 < t < 1100 s. Figure A4. Open in new tabDownload slide Same as for Fig. A2 for Pc = 60 MPa, except in (a) here the stopping criteria is the number of inversion iterations and the verification that the objective function becomes smaller than 10−3; (b) a profile is plotted each ∼49 s. Figure A4. Open in new tabDownload slide Same as for Fig. A2 for Pc = 60 MPa, except in (a) here the stopping criteria is the number of inversion iterations and the verification that the objective function becomes smaller than 10−3; (b) a profile is plotted each ∼49 s. Figure A5. Open in new tabDownload slide Same as for Fig. A2 for Pc = 95 MPa, except in (a) stopping criteria: 1500 iterations and the verification that the objective function becomes smaller than 10−3; (b) a profile is plotted each ∼175 s. Figure A5. Open in new tabDownload slide Same as for Fig. A2 for Pc = 95 MPa, except in (a) stopping criteria: 1500 iterations and the verification that the objective function becomes smaller than 10−3; (b) a profile is plotted each ∼175 s. APPENDIX B: DETAILED RESULTS OF THE DETERMINISTIC INVERSION We present in this section the objective function as well as the reconstructed model of pore pressure along a longitudinal section of the fault that intersects the two boreholes issued from the deterministic inversion for the different injection tests at 30, 60 and 95 MPa confining pressure (Figs A2, A4 and A5). In the reconstructed model, the pressure at the injection borehole is by definition imposed, and no fluid flows outside the fault boundaries (grey areas). Finally, Fig. A3 shows a zoom version of Fig. 3(a) over different time intervals. © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Fault’s hydraulic diffusivity enhancement during injection induced fault reactivation: application of pore pressure diffusion inversions to laboratory injection experiments
JF - Geophysical Journal International
DO - 10.1093/gji/ggaa446
DA - 2020-10-15
UR - https://www.deepdyve.com/lp/oxford-university-press/fault-s-hydraulic-diffusivity-enhancement-during-injection-induced-qab0nkqZvu
SP - 2117
EP - 2132
VL - 223
IS - 3
DP - DeepDyve
ER -