Modeling of Fiber Optic Strain Responses to Shear Deformation of Fractures

Abstract

Identifying distributed strain sensing (DSS) patterns (or signatures), particularly those arising from different hydraulic fracture geometries, has gained significant attention and research effort. Recent works have generated a catalogue of signatures for planar hydraulic fractures in an elastic rock formation. Yet, in numerous cases (e.g., fault motion and some geothermal reservoir stimulation), the main mode of deformation is a shear on a fracture or a network of natural fractures (particularly during low pressure injection/circulation). However, the specific fiber signatures that result from such shear deformation have not been studied. In this study, we use a three-dimensional poroelastic hydraulic fracture simulator to capture the strain signatures resulting from the shear deformation of fractures in various orientations with respect to the monitoring well. Five key cases are examined: one where the fracture strike is perpendicular to the fiber, another with the strike running parallel to the fiber, a third case where the fracture strike is at 45 degrees to the fiber, a fourth case with a strike slip fault perpendicular to the fiber, and a fifth case where fiber is intersecting the fracture. Theoretically and physically meaningful results were obtained in all five cases, which completely differ from the heart-shaped signature of tensile fracture propagation. It was discovered that the strain pattern changes with the shear deformation direction with respect to the fiber. The model is then used to simulate the response of a fracture network at Utah FORGE to injection to assess whether a signature might be expected in response to the planned injection and circulation rates, and, if so, what strain pattern might be expected. The simulation confirms that a strain response can indeed be observed. More importantly, the fiber response that would be detected in the monitoring well would be a combination of strain signatures from dilation and shear deformation of differently oriented natural fractures. The results in this study provide useful insights on the application of fiber to other stimulation and/or circulation scenarios where shear deformation of a fracture or fracture network plays a major role.

1. Introduction

Distributed strain sensing (DSS) technology has become an essential tool in unconventional petroleum reservoirs, and it has recently seen a surge in application in enhanced geothermal system (EGS) development. It has proven useful in refining estimates of stimulation outcomes based on traditional methods such as micro-seismic and pressure monitoring. In DSS operations, the linear strain experienced by the fiber at specific moments in time is measured [1,2]. By conducting repeated measurements at closely spaced time intervals, it is possible to construct a detailed strain signature. A standard DSS setup includes a laser interrogator that emits laser signals, captures the light that is backscattered by the geological formation, and translates these data into strain measurements at designated points along the fiber [3]. The spacing between these points, known as the gauge length, is a critical parameter [4]. In comparison with other downhole sensors like geophones, DSS provides information over a longer length and offers real-time data acquisition.

The DSS method measures linear strain variations across the fiber optic cable for a specified gauge length [5]. The calculation of strain (${\epsilon }_{xx}^i$) and strain rate (${\dot{\epsilon }}_x^i$) along the fiber uses fundamental principles of strain from continuum mechanics:

$${\epsilon }_{xx}^i=\frac{u_x^{i+1}-u_x^i}{L}$$

(1)

$${\dot{\epsilon }}_x^i=\frac{d{\epsilon }_x}{dt}=\frac{{\epsilon }_x^{n+1}-{\epsilon }_x^n}{t^{n+1}-t^n}$$

(2)

where i and n are the location of the sensing point and time (t) index, respectively, x is the direction along the fiber cable, u_x indicates displacement in the direction along the fiber, and L is the gauge length (Figure 1).

There is a keen interest in identifying the strain patterns generated from both single planar fractures and complex networks of fractures. Much effort has been devoted to developing a comprehensive catalogue of signatures corresponding to single and multiple propagating hydraulic fractures (e.g., [6,7,8,9,10]) using elastic analysis of planar fractures. Yet, there is a noticeable gap in the literature concerning the strain signatures with regard to shear deformation of fractures, which is of major interest in enhanced geothermal systems (EGS) and other geothermal applications. In this research, we use a 3D poroelastic fracture simulator that uses a coupled displacement discontinuity (DD) method with the finite element method for fluid flow. This fracture simulator is used to identify the strain patterns that are induced at the fiber by the shear deformation of fractures for geothermal as well as carbon storage and other subsurface storage needs.

2. Materials and Methods

To analyze the distribution of fiber strain, we use the 3D distribution of deformation related to fracture deformation using a three-dimensional poroelastic displacement discontinuity method [11,12,13] with finite element for fluid flow in the fracture [14,15,16]. The model’s main components are highlighted below.

The governing equations of poroelasticity are obtained by combining constitutive equations and conservation laws [17]. These governing partial differential equations can be solved using a displacement discontinuity method with a boundary collocation technique. The fracture is discretized into a number of elements (N) and the total time is divided into number of time steps (p + 1). The induced normal stress, shear stresses, and pore pressures at time t, on element “J”, due to a constant spatial distribution of continuous normal and shear DDs and sources on all other elements are given by [14,18,19,20,21]:

$$\begin{array}{c}{\sigma }_{i2}^J\left(\mathit{X},t\right)={\displaystyle \sum _{r=1}^N}\left[{\sigma }_{i2k2}^{cd}\left(\mathit{X}-\mathit{\chi },t-{\tau }_p\right)\mathsf{\Delta }{\psi }_{k2}^{rp}+{\sigma }_{i2}^{cs}\left(\mathit{X}-\mathit{\chi },t-{\tau }_p\right)\mathsf{\Delta }{\varphi }^{rp}\right]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{l=0}^{p-1}}\sum _{r=1}^N\left[{\sigma }_{i2k2}^{cd}\left(\mathit{X}-\mathit{\chi },t-{\tau }_l\right)\mathsf{\Delta }{\psi }_{k2}^{rl}+{\sigma }_{i2}^{cs}\left(\mathit{X}-\mathit{\chi },t-{\tau }_l\right)\mathsf{\Delta }{\varphi }^{rl}\right]\end{array}$$

(3)

$$\begin{array}{c}{\mathit{p}}^{\mathit{J}}\left(\mathit{X},\mathit{t}\right)={\displaystyle \sum _{\mathit{r}=1}^{\mathit{N}}}\left[{\mathit{p}}_{\mathit{k}2}^{\mathit{c}\mathit{d}}\left(\mathit{X}-\mathit{\chi },\mathit{t}-{\mathit{\tau }}_{\mathit{p}}\right)\mathsf{\Delta }{\mathit{\psi }}_{\mathit{k}2}^{\mathit{r}\mathit{p}}+{\mathit{p}}^{\mathit{c}\mathit{s}}\left(\mathit{X}-\mathit{\chi },\mathit{t}-{\mathit{\tau }}_{\mathit{p}}\right)\mathsf{\Delta }{\mathit{\varphi }}^{\mathit{r}\mathit{p}}\right]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{\mathit{l}=0}^{\mathit{p}-1}}{\displaystyle \sum _{\mathit{r}=1}^{\mathit{N}}}\left[{\mathit{p}}_{\mathit{k}2}^{\mathit{c}\mathit{d}}\left(\mathit{X}-\mathit{\chi },\mathit{t}-{\mathit{\tau }}_{\mathit{l}}\right)\mathsf{\Delta }{\mathit{\psi }}_{\mathit{k}2}^{\mathit{r}\mathit{l}}+{\mathit{p}}^{\mathit{c}\mathit{s}}\left(\mathit{X}-\mathit{\chi },\mathit{t}-{\mathit{\tau }}_{\mathit{l}}\right)\mathsf{\Delta }{\mathit{\varphi }}^{\mathit{r}\mathit{l}}\right]\end{array}$$

(4)

where p + 1 is the number of time increments used (time steps), N is the number of elements, $\mathsf{\Delta }{\psi }_{k2}^{rl}$ and $\mathsf{\Delta }{\varphi }^{rl}$ are the increments of strengths of DDs and sources that occur on element r, and subscript k denotes the direction and can be 1 or 2. ${\sigma }_{i2k2}^{cd}$ and $p_{k2}^{cd}$ are the boundary influence coefficients and they denote the actual stress and pore pressure at the midpoint of element M due to a constant DD applied on element N.

The above Equations (3) and (4) form a set of four linear equations with four unknowns of one normal DD, two shear DDs, and fluid source. These unknowns can be found if the normal stress, two shear stresses, and pressure on the fracture are prescribed. However, in a fluid injection procedure, these stresses and pressure on the fracture surface tend to change with time and are not known beforehand so the fracture deformation equations need be coupled to the fluid transport inside the fracture.

2.1. Fluid Flow in the Fracture

The model used in this study to simulate the fluid flow inside the fracture is the cubic law [22] in which it is assumed that the fracture consists of two parallel plates separated by the fracture aperture. Additionally, fluid (assumed incompressible) continuity law, which governs the fluid flow inside fractures, is also used. The following approximate solution (Equations (5)–(10)), based on the Galerkin finite element method, is determined by combining cubic law, fluid continuity equation and deriving the weak form of the resulting partial differential equation [14]:

$$\mathit{M}=\sum _{e=1}^N\left({\mathit{M}}_e\right), { \mathit{M}}_e={\int }_{-1}^{+1}{\int }_{-1}^{+1}{\mathit{N}}^{T}A\mathit{N}\mathit{J}d{\xi }_{1}d{\xi }_2$$

(5)

$$\mathit{K}=\sum _{e=1}^N\left({\mathit{K}}_e\right), { \mathit{K}}_e={\int }_{-1}^{+1}{\int }_{-1}^{+1}\mathit{B}\frac{w^3}{12\mu }{\mathit{B}}^T\mathit{J}d{\xi }_{1}d{\xi }_2$$

(6)

$$\mathit{F}=\sum _{e=1}^N\left({\mathit{F}}_e\right), { \mathit{F}}_e={\int }_{-1}^{+1}{\int }_{-1}^{+1}{\mathit{N}}^T\left(-2\nu \right)\mathit{J}d{\xi }_{1}d{\xi }_2$$

(7)

$$\mathit{B}=\mathit{D}\mathit{N}{\mathit{J}}^{-1}, \mathit{D}\mathit{N}=\partial N_i/\partial {\xi }_i$$

(8)

$$\mathit{Q}=\left[0.0\dots Q_{inj}\left(t\right)\dots -Q_{ext}\left(t\right)\dots 0.0\right]$$

(9)

$$\mathit{M}\dot{\mathit{p}}=\mathit{K}\left(\mathit{p}\right)\mathit{p}=\mathit{F}\left(\mathit{p}\right)$$

(10)

where N is the shape functions matrix, J is the Jacobian matrix, w is the fracture aperture, μ is the fluid viscosity, A is the area of each element, N is the number of elements, p is the pressure, $\dot{\mathit{p}}$ is the pressure derivative, K is the stiffness matrix, and Q is the matrix with injection or extraction rate at the nodes where injection or extraction happens.

Equations (5)–(10) relate to three unknowns: the pressure inside the fracture, fracture aperture, and leak-off from the fracture (denoted by F). Because the fracture aperture depends on the pressure inside the fracture, Equation (10) should be solved in an iterative manner along with DD (Equations (3) and (4)). First, the total time is divided into several time steps, and for the first time step, Equation (10) is solved assuming zero leak-off and initial fracture aperture, and then the calculated pressure is substituted in poroelastic DD (Equations (3) and (4)). Then, the amount of leak-off is updated and the new fracture aperture is calculated based on the normal and shear DD distribution. For the subsequent time steps, initial values of DDs and leak-off are obtained from the solution of the previous time step. This procedure is performed iteratively until the pressure inside the fracture converges within a prescribed error [14]. Combining DD equations and fluid flow equations, a unique solution can be obtained for all DD components (shear and normal), pressure, and stress distributions on fracture surfaces. Two shear DD components for each element is the shear deformation that particular element underwent in the given time step.

2.2. Modeling the Mechanical Behavior of Fractures

As indicated before, the DD method is used to model both natural (joint) and hydraulic fractures. The transition from one to another is dictated by the effective stress, i.e., if the natural fracture has zero effective stress (the total stress on the fracture walls is equal to the fluid pressure inside the fracture), it will be treated as a pressurized crack that is mechanically open. Moreover, a fracture that is open cannot support shear stress on its surfaces. If a fracture’s surfaces have not opened completely, it is treated as a joint and the surfaces can carry normal and shear stresses. In this condition, is assumed that the shear force on the fracture surface is a function of the shear DD in the same direction. For the joint fracture, the following constitutive equations are considered as the complementary to the discontinuity equations employed:

$${\mathit{\sigma }}_{\mathbf{33}}\left({\mathit{X}}^{\mathit{m}},{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}\right)\mathbf{=}{\mathit{f}}_{\mathbf{3}}\left({\mathit{D}\mathit{D}}_{\mathbf{33}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)\right)\mathbf{+}{\mathit{f}}_{\mathit{s}}\left({\mathit{D}\mathit{D}}_{\mathbf{13}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right),{\mathit{D}\mathit{D}}_{\mathbf{23}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)\right)$$

(11)

$${\mathit{\sigma }}_{\mathbf{13}}\left({\mathit{X}}^{\mathit{m}},{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}\right)\mathbf{=}{\mathit{f}}_{\mathbf{1}}\left({\mathit{D}\mathit{D}}_{\mathbf{13}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)\right)$$

(12)

$${\mathit{\sigma }}_{\mathbf{23}}\left({\mathit{X}}^{\mathit{m}},{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}\right)={\mathit{f}}_{\mathbf{2}}\left({\mathit{D}\mathit{D}}_{\mathbf{23}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)\right)$$

(13)

where ${\sigma }_{33}\left({\mathit{X}}^{\mathit{m}},t_{n+1}\right)$, ${\sigma }_{13}\left({\mathit{X}}^{\mathit{m}},t_{n+1}\right)$, and ${\sigma }_{23}\left({\mathit{X}}^{\mathit{m}},t_{n+1}\right)$ are the local normal, first local shear, and second local shear tractions on the center of the element (m) after $t_{n+1}$ time; $f_3\left({DD}_{33}^{t_{n+1}}\left({\mathit{X}}^{\mathit{m}}\right)\right)$ describes closing/opening of the joints; $f_s\left({DD}_{13}^{t_{n+1}}\left({\mathit{X}}^{\mathit{m}}\right),{DD}_{23}^{t_{n+1}}\left({\mathit{X}}^{\mathit{m}}\right)\right)$ explains the dilation effect; and $f_1\left({DD}_{13}^{t_{n+1}}\left({\mathit{X}}^{\mathit{m}}\right)\right)$, $f_2\left({DD}_{23}^{t_{n+1}}\left({\mathit{X}}^{\mathit{m}}\right)\right)$ define shear traction versus shear deformation relations. The normal behavior of a joint is solely related to the effective stress on the fracture surface. According to the fundamentals of joint-fracture behavior, the constitutive model in the normal direction can be represented as follows:

$${\mathit{\sigma }}_{\mathbf{33}}\left({\mathit{X}}^{\mathit{m}},{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}\right)\mathbf{=}{\mathit{\sigma }}_{\mathbf{33}}\left({\mathit{X}}^{\mathit{m}},{\mathit{t}}_{\mathit{n}}\right)\mathbf{+}\left(\mathsf{\Delta }{\mathit{\sigma }}_{\mathbf{33}}^{\mathit{M}\mathit{e}\mathit{c}\mathit{h}\mathit{a}\mathit{n}\mathit{i}\mathit{c}\mathit{a}\mathit{l}}\left({\mathit{X}}^{\mathit{m}},{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}\right)\mathbf{+}\mathsf{\Delta }{\mathit{\sigma }}_{\mathbf{33}}^{\mathit{D}\mathit{i}\mathit{l}\mathit{a}\mathit{t}\mathit{i}\mathit{o}\mathit{n}}\left({\mathit{X}}^{\mathit{m}},{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}\right)\right)$$

(14)

$$\begin{array}{c}\mathsf{\Delta }{\mathit{\sigma }}_{\mathbf{33}}^{\mathit{M}\mathit{e}\mathit{c}\mathit{h}\mathit{a}\mathit{n}\mathit{i}\mathit{c}\mathit{a}\mathit{l}}\left({\mathit{X}}^{\mathit{m}},{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}\right)\mathbf{+}\mathsf{\Delta }{\mathit{\sigma }}_{\mathbf{33}}^{\mathit{D}\mathit{i}\mathit{l}\mathit{a}\mathit{t}\mathit{i}\mathit{o}\mathit{n}}\left({\mathit{X}}^{\mathit{m}},{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\mathit{K}}_{\mathit{n}\mathit{o}\mathit{r}\mathit{m}\mathit{a}\mathit{l}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\times \mathsf{\Delta }{\mathit{D}\mathit{D}}_{\mathbf{33}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)-{\mathit{K}}_{\mathit{n}\mathit{o}\mathit{r}\mathit{m}\mathit{a}\mathit{l}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\times {\mathit{D}\mathit{D}}_{\mathit{d}\mathit{i}\mathit{l}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)\end{array}$$

(15)

$$\mathsf{\Delta }{\mathit{D}\mathit{D}}_{\mathbf{33}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)={\mathit{D}\mathit{D}}_{\mathbf{33}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)-{\mathit{D}\mathit{D}}_{\mathbf{33}}^{{\mathit{t}}_{\mathit{n}}}\left({\mathit{X}}^{\mathit{m}}\right)$$

(16)

$${\mathit{D}\mathit{D}}_{\mathit{d}\mathit{i}\mathit{l}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)=\mathsf{\Delta }{\mathit{D}\mathit{D}}_{\mathit{s}\mathit{h}\mathit{e}\mathit{a}\mathit{r}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)\times \mathit{t}\mathit{a}\mathit{n}{\mathit{\varphi }}_{\mathit{d}\mathit{i}\mathit{l}}$$

(17)

$${\mathit{D}\mathit{D}}_{\mathit{s}\mathit{h}\mathit{e}\mathit{a}\mathit{r}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)=\sqrt{{\left({\mathit{D}\mathit{D}}_{\mathbf{13}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)\right)}^{\mathbf{2}}\mathbf{+}{\left({\mathit{D}\mathit{D}}_{\mathbf{23}}^{{\mathit{t}}_{\mathit{n}\mathbf{+}\mathbf{1}}}\left({\mathit{X}}^{\mathit{m}}\right)\right)}^{\mathbf{2}}}$$

(18)

where ${\sigma }_{33}^{Mechanical}\left({\mathit{X}}^{\mathit{m}},t_{n+1}\right)$, and $\mathsf{\Delta }{\sigma }_{33}^{Dilation}\left({\mathit{X}}^{\mathit{m}},t_{n+1}\right)$ are the mechanical and dilation share of traction on the element (m), ${\mathit{K}}_{normal}^{t_{n+1}}$ is the normal stiffness of joint at time ${\mathit{t}}_{\mathit{n}+1}$, and ${\varphi }_{dil}$ is the dilation angle. The shear behavior of rock joints is modeled using elastic-plastic behavior (i.e., softening or hardening due to a decrease or increase in the effective stress, respectively). The Mohr–Coulomb failure criterion is used to distinguish the beginning of the plastic behavior. Because the maximum tolerable shear depends on the amount of effective stress, the choice of softening or hardening depends on the rate of change in the effective stress in the element. If the effective stress decreases, the element follows the softening branch during the poroelastic behavior. In contrast, if the effective stress increases but permanent slippage on the joint surface is present, the element follows the hardening branch [14].

2.3. Calculating Strain/Strain Rate

To calculate the strain and strain rate for each case, the HF simulator is executed first using the desired set of input data on fracturing conditions and the fiber location. At each step of propagation, reservoir displacements at fiber locations are obtained. Finally, strain and strain rates are calculated using Equations (1) and (2) and waterfall plots are obtained. A summary of this workflow is shown in Figure 1b.

3. Results

3.1. Strain Signatures in Response to Shear Deformation of a Fracture

In this paper, we consider five different ways that a natural fracture or a fault is positioned in relation to the fiber to show the strain patterns caused by shear deformation. The first scenario is for a fracture strike oriented perpendicular to the fiber and slipping in the direction of the fiber (dip slip), the second case is for a fracture strike parallel to the fiber and having dip slip, the third case considers a fracture strike at 45 degrees to the fiber and having dip slip, the fourth case is with a strike slip fault perpendicular to the fiber, and the fifth case is with fiber intersecting the fracture. The properties of the formation and other parameters used for this analysis are shown in

Table 1. Input data for cases 1–5 [23].

E	Young’s Modulus	5.387 × 104	MPa
ν	Drained Poisson’s ratio	0.29	-
νu	Undrained Poisson’s ratio	0.35	-
φ	Porosity	0.05	-
μ	Fluid viscosity	1.00 × 10−3	Pa.s
B	Skempton’s coefficient	0.47	-
cf	Fluid diffusivity	3.08 × 10−5	m2/s
α	Biot’s effective stress coefficient	0.69	-
k	Permeability	4.5 × 10−5	Darcy
φ	Friction angle	35	Degrees
c	Cohesion	0.6 × 106	Pa
Ks	Shear stiffness	2.0 × 1010	Pa/m
Kn	Normal stiffness	2.0 × 1010	Pa/m
Sv	Vertical stress	30.87	MPa
SH	Maximum horizontal stress	25.97	MPa
Sh	Minimum horizontal stress	24.89	MPa

. For all cases, two scenarios are considered: one that considers the strain signature due to only the effect of shear deformation, and the second that considers the effect of both shear deformation and dilation of the fracture. It should be noted that for cases 1–5, even though the effect of poroelasticity is taken into account via fracture deformation, the effect of pore pressure transience is not included in fiber strain calculation to allow a direct comparison with the expected signature of a tensile hydraulic fracture published in the literature.

3.1.1. Case 1: Natural Fracture Strike Is Perpendicular to the Fiber

In this scenario, we consider a rectangular-shaped fracture with the dimensions 20 m × 5.65 m. The fiber optic cable is placed 20 m away from the center of this fracture, as shown in Figure 2a. The fracture is discretized into a total of 560 elements, each with a size of 0.5 m × 0.6 m. The dip angle of the fracture is assumed to be 45 degrees. Next, fluid is injected into the fracture center at a small injection rate of 0.5 L/s. The injection is carried out for 60 min at a constant rate. A time step value of 1 min is used to record the strain in the fiber. The obtained results are presented in the subsequent figures.

Figure 3 shows the distribution of shear deformation (along dip) and dilation of the fracture at various times. It can be seen that as injection continues, the area of the fracture experiencing shear deformation increases.

Figure 4a displays the recorded strain signature along the fiber at different times as a result of the shear deformation of the fracture. It can be seen that the central region of the fiber undergoes tension while its ends experience compression (see also Figure 2b) and the results are clearly symmetrical around the 0 value of the y-axis. As a result of shear deformation, the rock mass on either side of the fracture moves in the opposite direction; because of this movement, the center of the fiber tends to stretch, and the two ends compress. The signature pattern of the combined effect of shear and dilation shown in Figure 4b is very similar to the shear only scenario but the magnitudes of strain are an order larger than the shear-only scenario. It is known that dilation of the fracture creates a tensile zone in the vicinity of the fracture tip (as does shear movement); therefore, the effect of dilation is the same as the effect of shear deformation for this fracture/fiber orientation, which explains the similarity of Figure 4a,b. Figure 5 shows the strain rate plots for only shear and dilation plus shear scenarios. The strain rate plot of shear plus dilation is similar to that of absolute strain, but this is not the case for the shear only case. According to the strain rate plot of Figure 5a, most of the tensile and compressional deformation in the fiber happens in the first 20 min of injection. This is because in the shear only case, the fracture slips in the first 20 min of injection and achieves a stable configuration with no further slip. In the dilation plus shear case, the fracture dilation which has the most effect on the fiber continues to increase as injection continues.

3.1.2. Case 2: Natural Fracture Strike Parallel to the Fiber

For this case, the fracture is 20 m × 5.65 m in size and the fiber extends for 60 m. The direction of the fracture dip is aligned in a way that shear deformation is towards the fiber. The fracture has a dip angle of 45 degrees. The positioning of both the fracture and the fiber are shown in Figure 6a. Fluid is injected into the fracture at a rate of 0.5 L/s and the resulting fiber strain response is calculated.

Figure 7 depicts how shear deformation and dilation change over time with continuous injection. It is noticeable that in the beginning, a section of the fracture near the injection well undergoes deformation. As the injection proceeds, the magnitude of this deformation increases. Also, the dilation area is significantly higher than shear deformation, which means dilation should have a higher effect on combined strain signature of shear and dilation. Figure 6b shows the shape of the fiber following the injection process. Upon analyzing the displacement magnitude at various points along the fiber, it is found that the fiber underwent deformation in both +x and −z directions. However, the deformation in the z-direction is minimal when compared to that of the x-direction, making it barely visible in the plot. According to the deformed shape, fiber is pushed away from the middle and therefore, it should be in tension due to the stretching. Figure 8a displays both the pattern and magnitudes of strain that were captured by the fiber as a result of shear deformation. As can be seen from the figure, the center portion of the fiber is in tension because of stretching due to the fracture shear deformation. When the middle portion of the fiber stretches, points located between −10, −20, and 10, 20 should undergo compression as the result of bending in the fiber. Figure 8b shows the combined signature of shear deformation and dilation. Similar to case 1, it can be seen that the two signatures are very similar in pattern, except for the fact that the magnitude of strains in shear and dilation scenario is significantly higher. Both dilation and shear deformation push the fiber away from the fracture and this combined effect should record a higher magnitude of strain.

Figure 8c,d show the strain rate plots corresponding to Figure 8a,b strain plots, respectively. These strain rate plots are quite similar to the plots for the absolute strain and the only significant difference is in the shear only plot in Figure 8c, where it is shown that most of the compressional strain at the two ends of the fiber occur during the first half of the injection. This is further confirmed in the shear deformation time plot shown in Figure 7a, where it is evident that most of the fracture slip takes place in the first 30 min of injection.

3.1.3. Case 3: Natural Fracture Strike Is 45 Degrees to the Fiber

In this example, the fracture dimensions are the same as those used in cases 1 and 2 (20 m × 5.65 m) with a dip angle of 45 degrees. However, the fracture is positioned in such a way that it strikes 45 degrees to the fiber. Shear deformation along the dip direction is considered in this case. Figure 9a provides a visual of the initial setup; the fracture’s center is at (0, 0, 0) and a 60 m long fiber is placed 20 m away from the fracture’s center. The fracture is injected at a flow rate of 0.5 L per second over a period of 60 min, and a one-minute time step is used for each measurement. The results of this case are presented below.

Figure 9b displays the deformed shape of the fiber and the direction of fracture deformation. As indicated by the arrows, the fracture shifts towards the fiber’s negative Y-direction and, as a result, the corresponding segment of the fiber is pushed away. Consistent with observations from earlier cases, the location of peak tensile deformation is where the fiber is bent, with compression ahead of it towards the Y = −30 m zone.

Figure 10a depicts the strain distribution from case 3, showing tension at the center and compression at both ends, as expected. Notably, the strain pattern in the negative Y section of the fiber is quite similar to that of case 2. Although there is some tension in the positive Y part of the fiber, its magnitude is comparatively small. Thus, it can be inferred that the side of the fiber that displays a greater compressive strain is in the direction of the fracture’s movement. This insight can be used to determine the direction of the potential shear deformation based on the observed strain signature. Figure 10b shows the combined strain signature of shear and dilation. It follows the same trend as cases 1 and 2 with the combined signature being similar to the shear only pattern but with a higher magnitude. The main reason for the similarity is that in all cases (1, 2 and 3), both shear deformation and dilation act in the same direction and their effects are aligned. Figure 11 shows the strain rate plots for shear only and shear plus dilation scenarios. Again, these plots are similar to that of strain plots in terms of the pattern. Further, the strain rate plot of shear only case shown in Figure 11a is not very smooth and the reason for this is the non-monotonic variability of the fracture shear deformation with time. This is why the irregular pattern is minimally visible in the shear plus dilation scenario where there is a continuous increase in fracture dilation, which has the most effect on strain signature.

For both parallel and perpendicular fractures to the fiber, simulations were repeated with increased injection rates to find out if the captured strain signature would be affected. It was found that such changes in the injection do not alter the overall pattern of the strain signature; rather, they increase the magnitude of the tensile and compressive regions observed. This increase in magnitude is attributed to the enhanced shear deformation resulting from a higher rate of fluid injection, which, in turn, leads to higher deformation of the fiber. Additionally, the impact of the fracture’s dip angle on the strain signature was explored, and it was found that as the dip angle increases, the magnitude of shear deformation increases, and as a result the magnitude of recorded strain increases. However, this does not affect the strain pattern. Therefore, the strain magnitude captured by the fiber is determined by the extent of shear deformation, while the strain pattern reflects the direction of the shear deformation and slip. This is very useful in identifying fault types and orientations from recorded fiber data.

3.1.4. Case 4: Strike Slip Fault Perpendicular to the Fiber

In this case, a strike slip fault perpendicular to the fiber is considered. The fault is sheared along the strike (X direction) and the fiber response is calculated. Fracture dimensions and injection rate are the same as the previous cases with injection carried out for 60 min. Figure 12a shows the initial configuration of the fiber and the fracture.

Figure 13 displays the fracture shear deformation and fracture aperture at selected times. Figure 14a shows the strain response recorded at the fiber for case 4 due to shear deformation. As can be seen from this figure, although the strain pattern is symmetrical around the Y = 0 axis, the strain sign is not. The largest compression of the fiber happens between Y values of 10, 30 and the largest tension occurs between Y values of −10, −30. A similar pattern can be observed between Y = 0, 10 and Y = 0, −10 but with a smaller magnitude. The reason for this behavior is that as the shear is applied and fracture starts to deform, the rock within the positive Y coordinate region is pushed towards the fiber while the rock in the negative Y coordinate zone is pulled away from the fiber. This behavior is clearly evident in Figure 12b.

Figure 14b shows the strain signature obtained for the combination of shear and dilation effects. In this case, shear deformation happens in a direction perpendicular to the fiber while dilation occurs parallel to the fiber, which means the two deformations act in different directions, leading to the differences observed in Figure 14a,b, unlike in cases 1, 2, and 3 where shear and dilation had the same effect on fiber. It is evident from the fracture aperture and shear deformation values shown in Figure 13 that the dilation effect is dominating. So, the shear deformation effect shown in Figure 14a is not manifested in the combined signature in Figure 14b. The strain rate plots for shear only and shear plus dilation cases are illustrated in Figure 13d and Figure 14c. The only difference between absolute strain and strain rate plots is the time it takes to show the deformation of fiber. This late response in strain plots (Figure 14a,b) is mainly due to the difference in magnitudes used to plot the contours.

3.1.5. Case 5: Fiber Intersects the Fracture

This case is similar to case 4 but instead of having the fiber at some distance from the fracture, it is placed across the fracture as illustrated in Figure 15a. The dimensions of the fracture and injection rates are consistent with the previous cases. The deformed shape of the fiber is shown in Figure 15b.

In all plots for case 5, a cutout can be seen from −1 to 1 m in the Y coordinates. This is because in the DD method, a numerical issue of singularity arises if the field points are too close to the fiber. Figure 16a shows the strain response recorded at the fiber for case 5 due to shear deformation. The pattern is slightly similar to that of case 4, but the compressional and tensile strain zones are closer to the fracture in case 5. Figure 16b shows the strain plot for the shear plus dilation case. As can be seen in this figure, the fiber records increasing compressional strain throughout the injection. The reason for this is as the fracture continues to open up with injection, the walls of the fracture push the fiber in opposite directions causing compression in the fiber. This phenomenon is similar to the compressional zone around a pressurized crack [23]. The magnitudes of strain continue to increase with time because the dilation of the feature keeps increasing with injection. Figure 16c,d show the strain rate plots that correspond to Figure 16a,b, respectively. It can be seen the strain rate is constant in both plots throughout the injection period, mainly due to the constant injection rate used.

3.2. Application to the FOGMORE Project of Utah FORGE

The main objective of this section is to assess the strain pattern resulting from a deforming fracture network during low pressure injection/circulation, as planned in the FOGMORE project (https://gdr.openei.org/submissions/1538; accessed on 12 September 2022) of Utah FORGE [24]. This is needed when studying reservoir stimulation involving DFN, as is the case for Well 16A in the Utah FORGE geothermal system, as well as when characterizing the network using low pressure circulation experiments. The main design question for the planned FOGMORE field experiment is whether the proposed injection rates can induce large enough rock mass deformation to be detected by a fiber in a monitoring well 100 m away. The Utah FORGE injection well 16A has had three stages of simulation and in this study, we focus on Stage 1, which is an open hole section that is 196 ft total length and intersects a number of natural fractures. The DFN used in the study is obtained from the Utah FORGE and is subject to revision based on the MEQ data.

It is expected that the pumping experiments during the planned circulation will not result in significant fracture propagation; rather, the pumping pressure will remain well below the S_h,min and thus, will only dilate the natural fracture with possible shear deformation. This process can be impacted by pore pressure diffusion and poroelastic effects, so we use a poroelastic displacement discontinuity method to calculate the strain in the monitoring production well in the Utah FORGE field. Figure 17 shows the position of the injection and production wells. The natural fractures that intersect the production well and form a network are also shown (using the DFN constructed in the Utah FORGE conceptual model of the reservoir). The production well is to be 100 m above the injection well based on current drilling plans. A fiber cable is to be placed in the production well as shown in Figure 17 and is used to record the strain caused by injection into Stage 1. The cable used for the simulation is assumed to be 55 m in length and the distance between two sensing points (i.e., the gauge length) is assumed to be 1 m. Therefore, a total of 56 sensing points are assumed to exist along the cable.

Table 2. Input data for cases 1–3 [23].

E	Young’s Modulus	5.387 × 104	MPa
ν	Drained Poisson’s ratio	0.29	-
νu	Undrained Poisson’s ratio	0.35	-
k	Permeability	4.5 × 10−5	Darcy
φ	Porosity	0.05	-
B	Skempton’s coefficient	0.47	-
μ	Fluid viscosity	1.00 × 10−3	Pa.s
cf	Fluid diffusivity	3.08 × 10−5	m2/s
α	Biot’s effective stress coefficient	0.69	-
Sv	Vertical stress	64.96	MPa
Sh	Minimum horizontal stress	45.11	MPa
SH	Maximum horizontal stress	55.77	MPa
-	Duration of injection	5	hours

shows the input properties used. Two scenarios are considered, using constant injection rates of 10 L/s and 1 L/s. The injection is carried out for a total of 5 h with a time step of 15 min.

Figure 18a,b show the calculated strain in the fiber with time due to the fracture network deformation because of fluid injection at rates of 10 L/s and 1 L/s into stimulation Stage 1. X is the distance measured along the fiber in the axis parallel to the production well, where x = 0 is the toe of the production well. As can be seen from the plot, recorded strain starts with a negative (tensile) value and becomes a positive (compressive) strain as injection continues. Figure 18c,d show the calculated strain rate for the two injection rates considered. It can be seen in both plots that there is a very short period of tensile strain in the fiber at the start of the injection and as injection continues, it records compression. There is a relatively high compression reading from 30 to 180 min and it gradually decreases towards the end of injection. Clearly, even a low injection rate leads to a strain signature so that the experiments can start from low rates and build to higher rates if necessary. It is interesting to notice that the strain/strain rate patterns do not match any patterns discussed in the single fracture cases considered above nor with patterns suggested by published catalogues. This can be explained by realizing the recorded strain is a combination of dilation and shear deformation of multiple natural fractures that are oriented in various directions. The dilation increases with time as the rock drains and becomes softer. As a result, while the entire dilating fracture remains closed, a portion experiences increasing ballooning which enhances the compressive stress shadow.

This can be better realized by considering the DFN aperture variations in time, as shown in Figure 19. The fractures with a high fracture opening contribute more by inducing a higher stress shadow in the surrounding rock. The apertures of fractures that are located above the injection well are higher than those of fractures below the well. This is because the fracture network portion below the well intersects with the wellbore at just one location while the upper part of the fracture network is intersected by the wellbore at multiple points and thus takes more fluid.

Figure 20 and Figure 21 show the shear deformation along the dip and strike of fractures in the x- and y- directions, respectively. As can be seen from these figures, there is some shear deformation on the fractures, but the magnitude is low when compared to that of fracture dilation. Because of these comparatively lower magnitudes, shear deformation’s contribution to the strain signature is minimal.

4. Conclusions

In this research, we utilized a 3D poroelastic displacement discontinuity (DD) model coupled with a finite element model for fluid flow to explore the potential strain patterns arising from the shear deformation of both a single fracture and a network of fractures. Five different scenarios, each involving various fracture/fault orientations with respect to the fiber, were considered, and strain catalogues were constructed. In case 1, where the fracture runs perpendicular to the fiber, shear deformation caused the two ends of the fiber to compress while the middle section experienced tension. In case 2, where the fracture is aligned parallel to the fiber, the trend was somewhat similar to case 1 but the strain pattern was different. The center of the fiber was in tension while the two ends were at compression, which was caused by the fiber being pushed away by shear deformation. In case 3, the fracture was oriented at a 45-degree angle to the fiber, and the resulting strain pattern in the negative Y portion of the fiber was similar to that of case 2. It was concluded that the side of the fiber that displays a similar pattern to case 2 is the direction that the fracture slip happens. A strike slip fault perpendicular to the fiber was considered as the fourth case, and, as expected, results showed that one portion of the fiber was being pulled while the other half was pushed away as a result of shear slip. In case 5, a scenario where the fiber intersects the fracture was considered and shear only results were quite similar to case 4, but shear plus dilation results were different because of the compression of the fiber due to the opening of fracture walls. It was also discovered that an increased injection rate increased the magnitude of recorded strain in all the cases; however, it did not affect the strain pattern. The model was then used to simulate the strain response to injection into a DFN located around Stage 1 of Well 16A in Utah FORGE. The results confirm that a strain signature can be obtained with the proposed low injection rates in a monitoring well 100 m away. It was found that the resulting strain pattern does not match existing catalogues and reflects a combination of patterns resulting from the dilation and shear deformation of natural fractures. The new strain catalogues related to shear deformation developed in this study can be used to help characterize geothermal reservoir stimulation, as well as unconventional reservoir stimulation that involves interactions between hydraulic and natural fractures.