# Modeling of fault activation and seismicity by injection directly into a fault zone associated with hydraulic fracturing of shale-gas reservoirs,

J. Rutqvist, Antonio Pio Rinaldi, F. Cappa, G. Moridis

► To cite this version:

J. Rutqvist, Antonio Pio Rinaldi, F. Cappa, G. Moridis. Modeling of fault activation and seismicity by injection directly into a fault zone associated with hydraulic fracturing of shale-gas reservoirs, Journal of Petroleum Science and Engineering, 2015, 127, pp.377-386. 10.1016/j.petrol.2015.01.019. hal-01231183

HAL Id: hal-01231183

https://hal.science/hal-01231183

Submitted on 16 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Modeling of Fault Activation and Seismicity by Injection
Directly into a Fault Zone associated with Hydraulic Fracturing
of Shale-gas Reservoirs

Jonny Rutqvist*, Antonio P. Rinaldi, Frédéric Cappa, and George J. Moridis

Lawrence Berkeley National Laboratory, Earth Sciences Division, Berkeley, CA 94 720, USA

* Corresponding author. Tel.: +1-510-486-5432, fax.: +1-510-486-5432
E-mail address: Jrutqvist@lb.gov (J. Rutqvist)

December 19, 2014

ABSTRACT

We conducted three-dimensional coupled fluid-flow and geomechanical modeling of fault activation and seismicity associated with hydraulic fracturing stimulation of a shale-gas reservoir. We simulated a case in which a horizontal injection well intersects a steeply dipping fault, with hydraulic fracturing channeled within the fault, during a 3-hour hydraulic fracturing stage. Consistent with field observations, the simulation results show that shale-gas hydraulic fracturing along faults does not likely induce seismic events that could be felt on the ground surface, but rather results in numerous small microseismic events, as well as aseismic deformations along with the fracture propagation. The calculated seismic moment magnitudes ranged from about -2.0 to 0.5, except for one case assuming a very brittle fault with low residual shear strength, for which the magnitude was 2.3, an event that would likely go unnoticed or might be barely felt by humans at its epicenter. The calculated moment magnitudes showed a dependency on injection depth and fault dip. We attribute such dependency to variation in shear stress on the fault plane and associated variation in stress drop upon reactivation. Our simulations showed that at the end of the 3-hour injection, the rupture zone associated with tensile and shear failure extended to a maximum radius of about 200 m from the injection well. The results of this modeling study for steeply dipping faults at 1000 to 2500 m depth is in agreement with earlier studies and field observations showing that it is very unlikely that activation of a fault by shale-gas hydraulic fracturing at great depth (thousands of meters) could cause felt seismicity or create a new flow path (through fault rupture) that could reach shallow groundwater resources.
# 1. Introduction

The rapid increase in North American shale-gas energy production has been made possible through new technology development, including extended-reach horizontal drilling and multistage hydraulic-fracture stimulation. But these new technologies have also raised concerns related to a range of local environmental problems (Arthur et al. 2008; Zoback et al. 2010)). One concern, investigated in this study, is whether shale-gas hydraulic fracturing could activate faults and thereby cause seismicity, opening up flow paths for upward fluid leakage and possible contamination of shallow potable groundwater resources (Arthur et al. 2008; Zoback et al. 2010; Davies et al., 2013; Rutqvist et al., 2013).

A first modeling study to investigate the potential consequences of fault reactivation during shale-gas hydraulic fracturing operations was presented in Rutqvist et al. (2013). Consistent with field observations, the study showed that a hydraulic fracturing operation to stimulate a deep shale-gas reservoir could only give rise to limited fault rupture, along with the possibility of (unfelt) microseismicity. In another study, Flewelling et al. (2013) used injection data and elastic fracture volume and length relationships to bound fracture-height data from 12,000 hydrofracturing stimulations conducted across North America. The hydraulic fracturing data showed that all microseismic events occurred less than 600 m above well perforation, although most were very much closer, and the farthest were usually associated with faults. These studies indicated that shale-gas hydraulic fracturing at great
depth (thousands of meters) could not create flow paths for leakage to reach shallow
2 groundwater resources.

3

Studies have also concluded that the likelihood of inducing felt seismicity during shale-gas
hydraulic fracturing operations, while not to be ruled out completely, is extremely small
(National Research Council, 2012; Davies et al., 2013). Indeed, after hundreds of thousands
of shale-gas fracturing stages conducted to date, only three examples of felt seismicity have
been documented (Davies et al., 2013). In Lancashire County, UK, two seismic events of
Richter scale magnitude *M*<sub>L</sub> = 2.3 and 1.5 were likely induced by direct injection into a
fault zone that had not been previously mapped (De Pater and Baisch, 2011). In another
case at the Eola Field of Garvin County, Oklahoma, in January 2011 (Holland, 2011), there
was a clear temporal correlation between the time of stimulation and the occurrence of 43
earthquakes that ranged in magnitude from *M*<sub>D</sub> = 1.0 to 2.8 (*M*<sub>D</sub> is the duration magnitude).
Finally, the third case of felt seismicity occurred at Etsho and Kiwigan fields in Horn River,
Canada, where 19 events between *M*<sub>L</sub> = 2 and 3 occurred having a clear temporal
correlation with the shale-gas operation; the largest (and felt) event, occurring in May 2011,
had a magnitude of *M*<sub>L</sub> = 3.8 (BC Oil and Gas Commission, 2012; Davies et al. 2013). Each
of these three cases of felt seismicity have been associated with reactivation of faults.

19

The biggest modeling uncertainty in the previous fault-activation modeling by Rutqvist et
al. (2013) was a 2D simplification of the full 3D field settings. In 2D plane-strain
simulations, it is difficult to estimate a representative injection rate, and some assumptions
have to be made about the shape of the rupture area (e.g., circular with diameter equal to
2D rupture length), which affects the calculated seismic magnitude. In this study we conduct, for the first time, a full 3D model simulation of fault activation associated with shale-gas fracturing. In such a 3D model simulation, the exact injection rate from the 3D field is a direct model input, and the seismic magnitude can be evaluated directly from the calculated rupture area and mean slip without the model uncertainties inherent in a 2D simplification. In this new 3D modeling study, we simulate the case in which a horizontal well intersects a subvertical fault, which then can be reactivated by injection directly into the fault. In addition, we investigate some issues not addressed in the previous 2D modeling in Rutqvist et al. (2013), including how the results correlate with fault and injection depth, fault dip, and fault frictional properties. We conclude with a discussion relating our modeling results to field observations and attempt to explain under which conditions a shale-gas fracturing stimulation could induce a felt seismic event.

## 2. Model Setup

We adopted the modeling approach that was applied in the previous 2D modeling study in Rutqvist et al. (2013). That is, we used the coupled multiphase fluid-flow and geomechanical simulator TOUGH-FLAC (Rutqvist, 2011) to model water-injection and fault responses, and we applied seismological theories to estimate the corresponding seismic magnitude. The fault was modeled as a discrete feature using finite thickness elements having anisotropic elasto-plastic properties. Shear failure was governed by a Mohr-Coulomb constitutive model with strain-softening frictional strength properties, consistent with a seismological slip-weakening fault model (Cappa and Rutqvist, 2011).
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

This allowed us to model sudden (seismic) slip events and to estimate their seismic magnitude. The adopted modeling approach has also been extensively applied for modeling fault activation associated with underground CO₂ injection (e.g. Cappa and Rutqvist, 2012; Mazzoldi et al., 2012; Rinaldi et al., 2014).

The model domain and the material properties are presented in Fig. 1 and Table 1, respectively. We model a full 3D-geological system (x, y, z: 2 km×10 km×2 km) generally tuned towards conditions that could be encountered in the Marcellus shale-gas play in the Northeastern U.S. This includes model input of *in situ* stress, fluid pressure, temperature, material properties, and injection rates. In a base-case simulation, we adopt conditions consistent with areas where the Marcellus shale is located at a depth of about 2000 m (6562 ft). The model is representative of the Marcellus shale-gas play with a 30 m thick gas-bearing shale, bounded at the top and bottom by other low-permeability formations (such as inorganic gray shale and limestone). This system is intersected by a steeply dipping fault, which in the base case has a dip of 80°. We simulate a case in which the horizontal injection well intersects the fault, and we directly inject the fluid volume related to a 3-hour hydraulic fracturing stage directly into the fault.

We set the initial conditions assuming linear pore pressure and temperature gradients (9.81 MPa/km and 25°C/km, respectively), with constant hydraulic boundary conditions (i.e., open to fluid flow), except for the planes x = 0 and y = 0 where a no-flow condition is applied (Fig. 1). Mechanical boundary conditions are null displacement at x = 0 and y = 0 planes, and constant stress elsewhere. The initial stress field is selected to represent the
conditions at the Marcellus shale play as detailed and justified in Rutqvist et al. (2013). We first set the vertical stress gradient (maximum principal stress) to 26,487 Pa/m, corresponding to an overburden density of about 2700 kg/m³. We then consider the minimum principal stress to be horizontal and oriented parallel to the horizontal well, which would lead to vertical hydro-fractures perpendicular to the well, but which in this case follow the weak planes of the fault. This does also correspond to a normal faulting stress field, in which the minimum horizontal stress (and minimum principal stress) is directed normal to the strike of the fault. We set the magnitude of the initial minimum horizontal stress corresponding to a horizontal-over-vertical stress ratio of R = σ<sub>h</sub>/σ<sub>v</sub> = 0.6. There are uncertainties in the horizontal-over-vertical stress ratio and, as highlighted by Rutqvist et al. (2013), this ratio has an impact on the magnitude of fault shear activation. However, several sources (e.g. Cipolla et al., 2010) indicate a fracture closure stress of about 0.7psi/ft and this corresponds to a horizontal-over-vertical stress ratio of R = σ<sub>h</sub>/σ<sub>v</sub> = 0.6. In this study, we keep the horizontal-over-vertical stress ratio fixed at R = σ<sub>h</sub>/σ<sub>v</sub> = 0.6, but vary the depth of the system, which also means a variation in stress magnitude at the depth of the injection. The magnitude of the intermediate stress, which in this case of a normal faulting stress regime would be oriented parallel to the fault strike does not affect the potential for shear failure along the fault.

19

Another important parameter in our analysis is the shear strength of the fault and how it evolves along with the reactivation. Here, we use the strain-softening Mohr-Coulomb model, in which the coefficient of friction and cohesion decreases with slip, i.e., once the

20

21

22
peak shear strength is achieved and the fault slips, the cohesion drops to zero and the coefficient of friction drops to a residual value. In the numerical model, this is simulated by reducing the coefficient of friction and cohesion from peak to residual values over a plastic shear strain of $10^{-3}$ (Cappa and Rutqvist, 2011). In the base case, we use a coefficient of friction of $\mu = 0.6$, with a residual value (after slip) equal to $\mu_R = 0.4$, whereas the cohesion drops to zero from an initial value of 1 MPa. A larger difference between the peak and residual friction values represents a more brittle behavior that is expected to lead to a larger shear-stress drop and seismic event. The selection of the frictional coefficient parameters are also discussed and justified in Rutqvist et al. (2013), acknowledging that this is one possible set of reasonable values of the frictional coefficient. The fault shear strength and how it weakens with slip is defined by a set of parameters that are varied in this study

Other fault properties as well as properties of the shale listed in Table 1 are equivalent to those used and justified in Rutqvist et al. (2013). In this study, we assume that the fault is nearly impermeable (hydraulically indistinguishable from the host rock), though the permeability and porosity can increase as a result of fracturing and shear. We consider a nearly impermeable fault a realistic assumption in this case. As pointed out by Flewelling et al. (2013), hydrocarbons cannot accumulate where there are permeable faults serving as pathways for buoyant oil and gas to leak upward. A relevant example of an impermeable fault in shale is a fault zone in Opalinus Clay exposed at the Mont Terri Rock Laboratory, Switzerland (Croisé et al., 2004). This zone is several meters thick, consists of intensively fractured rocks, has an inferred shear offset of 5 m, but is still hydraulically indistinguishable from the host rock, having an estimated permeability $k = 2 \times 10^{-20} \text{ m}^2$
(Croisé et al., 2004). Consistent with such a conceptual model, we considered 3.3 m thick fault where fracturing can take place along weak planes such as healed fractures along the fault plane, meaning that the tensile strength is low and taken as 1.7 MPa, based on the extension of the Mohr-Coulomb failure envelope into tension. Such a low tensile strength is reasonable, considering published shale strength data, including laboratory measurements on bedding planes and tectonic shears in Opalinus Clay (Bock, 2009; Young et al., 2010), and past laboratory experiments on Marcellus black-shale samples (Heard and Lin, 1986).

We consider porosity and permeability changes in the fault with tensile and shear rupture, according to a conceptual model of fracturing and shear along a fault zone rather than opening of a single hydraulic fracture. In this model, as soon as tensile or shear failure propagates within the fault plane, the porosity increases with plastic tensile strain $\epsilon_t^p$ and plastic shear strain a $\epsilon_s^p$ according to:

$$ \phi = \phi_0 + \phi_f = \phi_0 + \epsilon_t^p + \epsilon_s^p \times \tan \psi \quad (1) $$

where $\psi$ is the shear dilation angle. The approach of changing porosity with rupture is associated with the constitutive ubiquitous joint model, where plastic tensile and shear strain can be extracted and used for superimposing porosity changes by the rupture. The porosity increase due to the fracturing provides additional fluid storage that is important to consider in the fluid-flow simulation and significantly affects the fracturing propagation.
Similarly, the tensile and shear rupture also provides permeability that is superimposed on the initial fault permeability. The adopted permeability-change model described in Rutqvist et al. (2013), considers changes in equivalent fractured rock permeability along the fault as a function of plastic strain normal to the fault plane:

$$k = k_0 + k_f = k_0 + A(\varepsilon_n - \varepsilon'_n)^3 \quad (2)$$

where $k_0$ is the initial fault permeability, A is a constant, and $\varepsilon'_n$ is a threshold strain related required crack opening displacement for onset of permeability changes. Here, we used $\varepsilon'_n = 1 \times 10^{-4}$ and $A = 1 \times 10^{-5}$, meaning that permeability would increase to about $1 \times 10^{-14}$ to $1 \times 10^{-11}$ m² for a plastic strain normal to the fault on the order of $1 \times 10^{-3}$ to $1 \times 10^{-2}$. This is a very substantial permeability change from an initial permeability of $1 \times 10^{-19}$ m², one that provides rapid pressure diffusion along the fractured elements with the fracture propagation.

We simulated a typical hydraulic fracturing stimulation stage with injection at a rate of 200 kg/s (75 bbl/min) for 3 hours. According to Fig. 1, the horizontal section of the well is explicitly represented and connected to elements at the fault to simulate injection focused at a packed-off section of the well. We assumed that a fracture had already been initiated and extended within the fault to about 10 m from the well. Such fracture initiation would in the field be the result of an initial pressurization to formation breakdown pressure (e.g. Mayerhofer et al., 2011). We assume this initial fracture because the model is not refined enough around the well to properly simulate fracture initiation and formation breakdown.
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

In the field, the injection is typically shut down for a few minutes before beginning the main stimulation stage (e.g., Mayerhofer et al., 2011), which in our simulation is assumed to be 3 hours long.

During the stimulation, the fracture is extended within the fault plane, as illustrated in Fig. 1. We then evaluate magnitude and frequency of microseismic events as well as the total extent of the rupture zone. The seismic moment $M_0$, and the moment magnitude $M_w$ are evaluated for newly ruptured patches on the fault (related to the propagation of the fracturing front). According to Kanamori and Anderson (1975) and Kanamori and Brodsky (2001), we have $M_0 = A\mu d$ and $M_w = (\log_{10} M_0 / 1.5) - 6.1$, where A [m$^2$] is the area of the rupture, $\mu$ [Pa] is shear modulus of the rock, and d [m] is the average slip on the rupture. In the modeling, a total of 2160 m$^3$ (571,611 gallons) of water were injected during the 3-hour simulated injection.

### 3. Modeling results

In the following subsections, we first present the results of the base-case simulation, considering injection at 2000 m depth into a fault dipping 80 degrees, and having a coefficient of friction $\mu = 0.6$ that upon reactivation drops to a residual value $\mu_R = 0.4$. We then present parameter studies to evaluate the sensitivity of the results to changes in injection depth (1000 to 2500 m), fault dip (75° to 90°), and fault friction coefficients (0.2 to 0.8). In each case, we evaluate the results in terms of number and maximum magnitude of microseismic events, and the total vertical rupture length along the fault.
1

2

**3.1 Base-case simulation results**

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

Fig. 2 presents the simulation results for the base case: Fig. 2a shows the pressure evolution, whereas Fig. 2b shows the magnitude of fault slip (the along-dip offset displacement between the hanging wall and foot wall of the fault) and moment magnitudes of seismic events. Note that the initial stress normal to the fault at the depth of the injection is about 32.5 MPa. Fig. 2a shows how the injection pressure increases to reach a peak value of about 43 MPa and then slowly decreases along with the failure propagation. Fig. 2b shows that shear displacement is initiated after a few minutes of injection—when the pressure exceeds the stress normal to the fault—and then slips gradually along with the propagation of the failure zone. The biggest microseismic events occur during the first few minutes of injection, and then microseismicity continues during the 3-hour injection at reduced magnitudes and frequency—but also continues for at least 10 minutes after terminating the injection (after 180 minutes). The simulated magnitude of the microseismic events generally ranges between $M_w = -2$ and 0. Fig. 2c to e show that at the end of the 3-hour injection, the failure zone has extended to a radius of almost 200 m from the well upwards and laterally, while only about 100 m downwards. This can be explained by the vertical gradient in stress that makes it easier to propagate a fracture upwards. Moreover, at the end of the 3-hour injection, the fault slip is almost 5 cm at the center of the failure zone (near the injection well). The simulation in Fig. 2 captures results typical for hydraulic fracturing, including (1) small pressure fluctuations, especially during the first 20 minutes when failure and shear slip occur relatively close to the injection well, (2) a gradual reduction in the pressure required to propagate the fracture for a given injection rate, and
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
finally, (3) a typical pressure fall-off after terminating the injection. After terminating the injection and shutting in the well (at 180 min), the well pressure declines to below the stress normal to the fault plane (32.5 MPa), in a shut-in pressure behavior that could be used in the field to estimate the minimum principal stress magnitude.

Overall, the result from the base-case simulation shows that maximum moment magnitude is less than 0.5, and the total rupture is less than 200 m up from the injection. Field observations show that during typical stimulation operations in the North American shale-gas basins (Marcellus, Barnett, Eagle Ford, Woodford, Haynesville, Horn River), the maximum moment magnitude in a stimulation stage ranges between $M_W$ -2.5 to (in rare cases) up to $M_W$ = 1.0, with the high-end values usually associated with faults (Warpinski, 2009; Warpinski et al., 2012). The simulated microseismic events are within this range, with the lower limit probably affected by the minimum size of element patches in the fault (10 by 10 m) that can rupture in one instance. This means that in our modeling, we model the relatively larger microseismic events occurring from shear slip along the fault plane, whereas in the field, there are numerous smaller-magnitude events perhaps occurring as a result of slip in small-scale fractures in the host rock surrounding the fault, and these are not resolved in our modeling.

## 3.2 Sensitivity to injection depth

In their review of thousands of hydraulic fracturing treatments in North American shale basins, Warpinski et al. (2012) noticed an increased microseismic magnitude and frequency with depth. In Fig. 3, we present an analysis in which we varied the injection depth from
1000 to 2500 m (3281 to 8202 ft). The trend is clear and shows an increasing maximum magnitude and frequency with depth (Fig. 3f). Note that, from Fig. 3a, b, and d, although the total shear slip magnitude increases with depth, the extent of the rupture zone is equivalent, i.e., always about 200 m up from the injection point. We attribute the increased total shear-slip magnitude with depth to a higher shear stress on the fault when going deeper. A higher shear results in a greater shear-stress drop upon reactivation, and the shear slip is proportional to the shear-stress drop for a given fracture extent and shear modulus. The higher shear stress can also explain the increased moment magnitude of seismic events with depth. Overall, such depth dependency is consistent with field observations in Warpinski et al. (2012).

11

## 3.3 Sensitivity to fault dip

Fig. 4 presents the results for a different fault dip, between 75 and 90°. The results in Fig. 4f show that the maximum magnitude and frequency increase for decreasing dip angle. The results in Fig. 4 also show a significant difference in tensile and shear behavior. The plastic tensile strain in Figure 4a (signifying tensile failure) propagates farthest for a vertical fault (90°), whereas the plastic shear strain in Fig. 4b (signifying shear slip) propagates farthest in the case of inclined faults (e.g., 75°). The tensile failure propagates the longest for a vertical fault because the initial stress normal to the fault is the lowest in this case (32.4 MPa), slightly lower than for the inclined faults (e.g., 33.2 MPa for the 75° fault). The very small events ($M_w = -3$ to $-2$) for a vertical fracture resulted from minute shear displacement (about 1 micron) when new surface patches failed in tension. Shear failure propagates
longer for inclined faults, as well as both shear slip and microseismic magnitudes increase,
because the initial shear stress is much higher (e.g., 5.3 MPa for the 75° fault compared to 0
for the 90° fault). However, even for an inclined fault, when the shear stress is the highest,
the event magnitudes are less than about M<sub>w</sub> = 0.7.

3.4 Sensitivity to slip-weakening model parameters

The coefficient of friction and its evolution during failure may affect the amount of stress
drop and thereby the amount of slip and seismic magnitude. In the base case, coefficient of
frictions μ = 0.6 and μ<sub>R</sub> = 0.4 were used. A coefficient of friction of 0.6 is commonly
applied in fault-stability analyses and has been identified as a lower limit value for the most
common rocks (Zoback, 2007). However, clay-rich fault rock and fractures in shale could
have a much lower coefficient of friction, especially under wet conditions (Zoback, 2007;
Samuelson et al., 2012). Another data source are the investigations of Opalinus Clay at the
Mont Terri underground research laboratory, which indicates a coefficient of friction of
about 0.4 along bedding planes (Bock, 2009; Young et al., 2010). One additional
complicating factor is that the coefficient of friction is rate dependent and is also dependent
of the clay content (Kohli and Zoback, 2013). Kohli and Zoback (2013) determined
coefficient of friction on fractures in Barnett, Haynesville, and Eagleford shale samples to
be in the range of 0.4 to 0.8, for high to low clay content. Moreover, unstable slip behavior
was noted for shale with lower clay content, i.e. when the initial coefficient of friction was
in the range of 0.6 to 0.8.

22
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

Because of the uncertainty in the selection of coefficient of friction, we conducted sensitivity studies involving both peak and residual friction values. Moreover, we varied the rate of frictional strength drop with plastic shear strain. In the base case it was assumed that the coefficient of friction drops from $\mu = 0.6$ to $\mu_R = 0.4$ over a plastic shear strain of $10^{-3}$.

As discussed in Cappa and Rutqvist (2011), shear strength drop with plastic shear strain can be related to a fault slip weakening model, commonly applied in earthquake seismology, and the evolution of shear stress with slip (Kanamori and Bodsky, 2004). A plastic shear strain of $10^{-3}$ across the fault width of 3.33 m corresponds to a fault shear slip of $3.33 \times 10^{-3}$ m (3.33 mm). This would correspond to the critical slip distance, $D_c = 3.33 \times 10^{-3}$ m in a slip weakening model, meaning that the coefficient of friction reduces to its residual value (dynamic friction) over a slip displacement of 3.33 mm (Aochi et al., 2014). In the base case, the maximum shear strain after the 3-hour injection is on the order of $10^{-2}$ and a shear maximum slip is up to 5 cm (Fig. 2). However, the incremental plastic shear strain and slip during each event is much smaller, with slip magnitudes on the order of $1 \times 10^{-5}$ to $5 \times 10^{-4}$ m for event magnitudes ranging from about $M_w = -2$ to 0. This means that the full frictional strength drop from $\mu = 0.6$ to $\mu_R = 0.4$ will not occur and this may in turn affect the shear stress drop during an event. Although a critical slip distance of 3.3 mm is within the range of laboratory measurements on rough fractures (Ohnaka, 2013), the value could vary by orders of magnitude depending on scale (Ohnaka 2003) and might be strongly heterogonous in a fault plane (Aochi et al., 2014). Moreover, recent laboratory data on shale by Kohli and Zoback (2013) indicated values on the order of tens of microns, i.e. very
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

small values. Again, the values could vary by orders of magnitude when going from core
sample scale to field scale (Ohnaka, 2003), leading to considerable uncertainty in this
parameter.

Acknowledging the uncertainties we investigate the sensitivity to changes in peak and
residual friction coefficient and the critical slip distance. Figure 5f shows that there is a
trend of fewer but higher magnitude microseismic events if reducing the critical slip
distance to as low as 10⁻⁵ m. With such a low critical slip distance, the full shear strength
drop corresponding to a drop in the coefficient of friction from μ = 0.6 to μR = 0.4 would
occur as soon as shear failure occurs. Moreover, quite surprisingly, the analysis showed the
rather small influence of peak and residual friction coefficient, for variations in the initial
coefficient of friction between μ = 0.5 to 0.8 and for residual coefficient of friction between
μR = 0.3 to 0.6. It is only when the residual coefficient of friction is reduced to 0.2 that we
see more significant change in calculated moment magnitudes and rupture length, including
a maximum moment magnitude of M<sub>W</sub> = 2.3 (Figure 6f). These relatively larger events can
be observed in Fig. 5d as more significant slip displacements occur at 70, 100, and 160
minutes. The M<sub>W</sub> = 2.3 event is of similar magnitude as those felt at the aforementioned
U.K. and Oklahoma incidents, though such small magnitude events may likely go
unnoticed. Note that if the residual coefficient of friction is μR = 0.2, then the shear stress
on the fault from the tectonic stress field (about 3.6 MPa) is higher than the residual shear
strength of the fault (about 2.5 MPa). This means that once fault activation is initiated,
shear strength within the rupture area will drop to below the prevailing shear stress, and
because of this weakening, the shear slip can self-propagate outside the zone of
1 pressurization. This is the reason why in the case of $\mu_R = 0.2$, the shear rupture is somewhat
2 larger. Nevertheless, this result indicates that the only way to produce a microseismic event
3 above $M_W = 2$ in the current model setting is for the residual fault shear strength to be less
4 than the prevailing shear stress.

5

## 4. Discussion

6

7
8 Overall, the results of this 3D simulation study are in agreement with the findings in the
9 previous 2D study in Rutqvist et al. (2013), in terms of potential extent of the rupture and
10 magnitude of seismic events. The moment magnitudes calculated in this study ranged from
11 about $M_W = -2.5$ to 0.5, except for the special case of very low residual coefficient of
12 friction of 0.2, when a magnitude $M_W = 2.3$ event was calculated. Over the duration of the
13 3-hour injection, repeated events and aseismic slip were shown to amount to up to 0.06 m,
14 with the total radius of rupture extending up to 200 m.

15

16 Consistent with field monitoring, the modeling shows that when faults are present,
17 somewhat larger seismic events are possible—compared to those associated with regular
18 hydraulic fracturing seismic events—because larger surface areas are available for
19 rupturing (Fisher and Warpinski, 2011: Warpinski et al., 2012). Indeed, regular hydraulic
20 fracturing seismic events average around $M_W = -2.5$, whereas higher than usual events up to
21 about $M_W = 0.5$ are usually associated with fracturing along faults (Warpinski et al., 2012).
22 Moreover, the modeling results of repeated small microseismic events and aseismic slip
23 and fracture opening is consistent with field observations—that the energy emitted from
microseismic events represents only a small fraction of the energy input or the energy to open the fracture (or fracture network) (Warpinski et al., 2012). Natural fractures present in shale (Gale et al., 2014) may also open aseismically by increased fluid pressure and associated reduction of effective stress (Rutqvist, 2014). In this context, weakly cemented fractures being common in shales (Gale et al., 2014) could have a permeability that is extremely stress dependent. Also consistent with our modeling, field observations of fracturing along faults have shown that microseismic events can travel several hundred meters upwards (Fisher and Warpinski, 2012; Flewelling et al., 2013; Yang and Zoback, 2014).

Our analysis showed that the frequency and magnitude of the microseismic events are dependent on the initial shear stress on the fault plane, while the total rupture length is limited by the volume of water injected. In the case of injection directly into the fault plane, the shear activation (shear failure) occurs simultaneously with the hydraulic fracturing (tensile failure), and these events takes place starting near the well and propagate away from the well bore in repeated microseismic events. At greater depth and for inclined faults that are more optimally oriented for higher shear stress, these microseismic events are of slightly higher magnitude. This depth dependency is also consistent with field observations at North American shale gas basins (Warpinski et al., 2012). An explanation from the modeling is that upon reactivation, shear stress is relieved, and this shear-stress drop is greater in the case of a higher initial shear stress, which in turn can results in larger slip for a given rupture area.
The total rupture length obtained in this study (about 200 m radius), will depend on the total injection volume and leak-off from the main hydraulic fracturing into surrounding rock, as well as the initial gas saturation. In this case the conceptual model is a 30 m thick gas bearing formation surrounded by gray shale and limestone layers that does not contain gas, and in this analysis we simplified the system to be fully water saturated. If considering some gas saturation in the gas bearing formation (e.g. 50%), the hydraulic fracturing process could be affected by the fact that gas is much more compressible than water therefore delaying the pressurization for a given injection rate (e.g. Rinaldi and Rutqvist, 2014). However, as shown in Kim et al. (2014), hydraulic fracturing of shale gas reservoirs may be affected by complex two-phase flow processes, including vertical gravity segregation within the created fractures. To investigate these effects Kim et al., (2014) conducted a hydraulic fracturing simulation of a shale-gas reservoir for the extreme cases of water saturated and gas saturated rock, and achieved similar radial extent of the stimulated fracture. Thus, we would not expect any major extension of the calculated rupture zone if considering some initial gas saturation in 30 m thick gas bearing formation. Nevertheless, the study in Kim et al., (2014) shows that the created fracture volume may not be equal the injected volume due to two-phase flow effects and fluid leak-off to the surrounding rock. Here, we injected 2160 m³ (571,611 gallons) of water over the 3-hour injection, which can be considered a typical average injection volume per fracturing stage (e.g., compared with data shown in Mayerhofer et al. (2011) and BC Oil and Gas Commission (2012)). If injection were to continue or were at a higher rate, the rupture would propagate slightly longer. However, this also means that the total length of the
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

rupture zone can be controlled by limiting the total injection volume and guided by seismic
monitoring.

One important question is, under what conditions could it be possible to induce felt seismic
events? We would emphasize that only three such cases have been documented among
hundreds of thousands of hydraulic fracturing treatments, meaning that these are very rare
events, and hence it might be difficult to produce such events even in a numerical model.

Our analysis for the conditions considered in this study indicates that the only way to
produce events on the order of magnitude M<sub>W</sub> = 2 to 3, is a very brittle fault with a residual
coefficient of friction lower than the prevailing shear stress on the fault (see results for a
residual coefficient μ<sub>R</sub> = 0.2 in Fig. 5). In this case, sufficient new surface area can be
ruptured in one instance to produce a seismic event on the order of M<sub>W</sub> = 2 to 3. Moreover,
consistent with field observations at the three known cases of felt events, these simulated
events occurred hours into the hydraulic fracturing stage, i.e., after 70, 100, and 160
minutes, in Fig. 5d. Each simulated event in Fig. 5d resulted in up to 1 cm fault slip, and the
total slip is greater than 6 cm, which would certainly put some strain on the horizontal well.

Deformations of the horizontal wells were also observed at the incidents at Lancashire
County, U.K. and Horn River, Canada. However, some of the shale-gas-fracturing-related
events reported in the literature, e.g., the events at the Eola Field, Oklahoma, occurred
many hours after shut-in of the injection and farther away from the injection wells. This
may indicate channeling of injection fluid that, after some time, then enters a critically
stressed fault that then can be reactivated. Such a scenario is more in line with the
conceptual model in previous 2D modeling in Rutqvist et al. (2013), where we considered a
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

hydraulic fracture propagating to intersect a fault, not injection directly into the fault. In that study, such a fault could be reactivated many hours after continuous injection to produce larger than usual microseismic events (e.g., magnitude 2). That is, if that fault had some permeability, were near critically stressed, and had brittle slip properties with a low residual strength. Thus, it could be that if fluid pressure migrates out of the shale layer it might encounter more brittle, faulted rock types that could be a source of larger than usual events.

## 5. Conclusions

We have conducted 3D model simulations of fault activation during a shale-gas hydraulic fracturing operation that eliminates some important model uncertainties associated with previous 2D modeling studies. We simulated a case tuned toward the conditions at the Marcellus shale simulation, a hydraulic fracturing stage of 2000 m, but with direct injection into a steeply dipping fault. Thus, the entire water volume from one hydraulic fracturing stage was injected to propagate a fracture along the fault. The results of the simulations are consistent with field observations regarding the possible microseismic magnitude and rupture length, indicating the occurrence of small (unfelt) microseismic events and fault ruptures that could propagate upwards a few hundred meters from the injection interval. The following main results were achieved in this modeling study:

21

22

23

* During the stimulation injection, the shear activation (shear failure) occurred simultaneously with the hydraulic fracturing (tensile failure), starting near the well and propagating away from the well in repeated microseismic events.
1. The moment magnitudes ranged from about $M_W$ = -2.5 to 0.5, except for the special case of a very brittle fault with a residual shear strength below the prevailing shear stress, when a magnitude $M_W$ = 2.3 event was calculated.

4. The microseismic magnitude increased with depth and for inclined faults that are more optimally oriented for higher shear stress, indicating the dependency of shear-stress drop on the seismic magnitude.

7. Over the course of the 3-hour injection, repeated events and aseismic slip amounted to up to 0.06 m, with the total radius of the shear rupture extending up to 200 m.

9. Overall, this 3D modeling study in agreement with previous findings showing that it is very unlikely that activation of a fault by shale-gas hydraulic fracturing at great depth (thousands of meters) could cause felt seismicity or create a new flow path (through fault rupture) that could reach shallow groundwater resources.

14. The analysis indicates that the very rare observations of felt seismicity associated with shale-gas fracturing (three documented out of several hundred-thousand fracturing stages) might have been caused by locally very unfavorable conditions, wherein fluid-pressure changes could reach a large section of a fault that was close to critically stressed for shear and having very brittle slip weakening properties in which the residual strength was lower than the prevailing shear stress.

20. Though our results in terms of seismic magnitudes are in agreement with field observations at various U.S. shale-gas plays exposed to various stress regimes, our analysis was for the
1

2

3

4

5

6

7

8

9

10

11

12

13

case of steeply dipping faults under a normal-faulting stress regime at 1000 to 2500 m depth. It might be different when operating a lot deeper such as could be the case in future shale gas production in China. In such case, the stress field will be much higher, the shale even tighter and potentially more brittle, and therefore potentially leading to larger magnitude events.

Nevertheless, adequate site characterization for identifying and avoiding faults should be a priority in any shale gas development. Moreover, continuous monitoring of induced seismicity from the start of the injection can be used to detect any runaway fracturing along faults. That is, to detect fracturing that propagates far above and below the gas-bearing shale layer, indicating potential reactivation of a fault.

## Acknowledgments

14

15

16

17

18

19

We thank the U.S. Environmental Protection Agency (EPA), Office of Research and Development, for supporting this study under an Interagency Agreement (DW-89-92378101) with the US Department of Energy at the Lawrence Berkeley National Laboratory. The views in this article are those of the authors and do not necessarily reflect the views or policies of the U.S. EPA. Mention of commercial software does not constitute endorsement by EPA or LBNL.

20

21

## References

22

23

Aochi H., Poisson B., Toussaint R., Rachez X., Schmittbuhl J., 2014. Self-induced seismicity due to fluid circulation along faults. Gephys J Int, 196, 1544–1563.
1

Arthur, J.D., Bohm B., Layne M., 2008. Hydraulic fracturing considerations for natural gas

2

wells of the Marcellus shale. Presented at the Ground Water Protection Council 2008

3

Annual Forum, Cincinnati, Ohio, USA, 21-24, Sep.

4

BC Oil and Gas Commission, 2012. Investigation of Observed Seismicity in the Horn River

5

Basin. <http://www.bcogc.ca/publications/reports.aspx>

6

Bock HR. 2009. A experiment: updated review of the rock mechanics properties of the

7

Opalinus Clay of the Mont Terri URL based on laboratory and field testing. Mont Terri

8

Project, Technical Report 2008-04.

9

Cappa F., and Rutqvist J., 2011. Impact of CO₂ geological sequestration on the nucleation

10

of earthquakes. Geophys Res Lett, 38, L17313, doi:10.1029/2011GL048487.

11

Cappa, F., Rutqvist, J., 2012. Seismic rupture and ground accelerations induced by CO₂

12

injection in the shallow crust. Geophysical Journal International, 190, 1784–1789.

13

Cipolla, C. L., Lolon, E.P., Erdle, J.C., Rubin, B., 2010. Reservoir modeling in shale-gas

14

reservoirs. Society of Petroleum Engineers (SPE) Reservoir Evaluation & Engineering.

15

August 2010, SPE 125530.

16

17

Croisé, J., Schlichenrieder, L., Marschall, P., Boisson, J.Y., Vogel, P., Yamamoto, S., 2004.

18

Hydrogeological investigations in a low permeability claystone formation: the Mont Terri

19

Rock Laboratory. Physics and Chemistry of the Earth 29: 3–15.

20

21

Davies R., Foulger G., Bindley A., Styles P., 2013. Induced seismicity and hydraulic

22

fracturing for the recovery of hydrocarbons. Marine and Petroleum Geology 45, 171-185.

23

De Pater, C.J., Baisch, S., 2011. Geomechanical study of Bowland shale seismicity.

24

Synthesis Report, November 2011. Cuadrilla Resources Ltd, UK.

25

26

Ferrill, D.A., McGinnis, R.N., Morris, A.P., Smart, K.J., Sickmann, Z.T., Bentz, M.,

27

Lehrmann, D., Evans M.A. 2014. Control of mechanical stratigraphy on bed-restricted

28

jointing and normal faulting: Eagle Ford Formation, south-central Texas. AAPG Bulletin,

29

v. 98, no. 11, p. 2165-2216.

30

31

Fisher, K., Warpinski, N., 2011. Hydraulic Fracture-Height Growth: Real Data. Paper was

32

prepared for presentation at the SPE Annual Technical Conference and Exhibition held in

33

Denver, Colorado, USA, 30 October–2 November 2011. Society of Petroleum Engineers,

34

SPE 145949.

35

36

Flewelling S.A., Tymchak M.P., Warpinski N. 2013. Hydraulic fracture height limits and

37

fault interactions in tight oil and gas formations. Geophysical Research Letter, VOL. 40,

38

3602–3606.

39
1 Gale, J. F. W., Laubach, S. E., Olson, J. E., Eichhubl, P., and Fall, A., 2014, Natural

2 fractures in shale: a review and new observations. AAPG Bulletin, v. 98, no. 11, p. 2165-

3 2216.

4

5 Holland, A., 2011. Examination of possibly induced seismicity from hydraulic fracturing in

6 the Eola Field, Garvin County, Oklahoma. Oklahoma Geological Survey Open-File Report

7 OF1-2011.

8

9 Kanamori, H., Anderson, D.L., 1975. Theoretical basis of some empirical relations in

10 seismology, Bull. Seism. Soc. Am., 65, 1073-1095.

11

12 Kanamori, H., Brodsky, E.E., 2001. The physics of earthquakes. Physics Today (June

13 2001), 34–40.

14

15 Kim J., Um E.S., Moridis G.J., 2014. Fracture propagation, fluid flow, and geomechanics

16 of water-based hydraulic fracturing in shale gas systems and electromagnetic geophysical

17 monitoring of fluid migration. Paper prepared for presentation at the SPE Hydraulic

18 Fracturing Technology Conference held in The Woodlands, Texas, USA, 4–6 February

19 2014. Society of Petroleum Engineers, SPE 168578. Kohli, A.H., Zoback, M.D. 2013.

20 Frictional properties of shale reservoir rocks, J. Geophys. Res. Solid Earth. 118, 1-17.

21

22 Mazzoldi, A., Rinaldi, A.P., Borgia, A., Rutqvist, J., 2012. Induced seismicity within

23 geologic carbon sequestration projects: Maximum earthquake magnitude and leakage

24 potential. International Journal of Greenhouse Gas Control, 10, 434–442.

25

26 Mayerhofer, M.J., Stegent, N.A. Barth J.O. and Ryan K.M. 2011. Integrating fracture

27 diagnostics and engineering data in the Marcellus Shale. presentation at the SPE Annual

28 Technical Conference and Exhibition held in Denver, Colorado, USA, 30 October–2

29 November 2011. SPE 145463.

30

31 National Research Council, 2012. Induced seismicity potential in energy technologies,

32 National Academies Press, Washington D.C., pp. 300.

33

34 Ohnaka M., 2003. A constitutive scaling law and a unified comprehension for frictional slip

35 failure, shear fracture of intact rock, and earthquake rupture. J Geophys Res, 108(B2),

36 2080.

37

38 Rinaldi A.P., Rutqvist J., and Cappa F. 2014. Geomechanical effects on CO2 leakage

39 through fault zones during large-scale underground injection. International Journal of

40 Greenhouse Gas Control. 20, 117–131.

41

42 Rinaldi A.P., Rutqvist J., Sonnenthal E., L., and Cladouhos T.T. 2014. Coupled THM

43 modeling of hydroshearing stimulation in tight fractured volcanic rock. Transport in Porous

44 Media (Published online April 2014). DOI 10.1007/s11242-014-0296-5.

45
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

Rutqvist, J., 2011. Status of the TOUGH-FLAC simulator and recent applications related to coupled fluid flow and crustal deformations. Computers & Geosciences, 37, 739–750.

Rutqvist J. Fractured rock stress-permeability relationships from in situ data and effects of temperature and chemical-mechanical couplings. Geofluids (Published online, 2014). doi: 10.1111/gfl.12089.

Rutqvist, J., Birkholzer, J.T., Cappa, F., Tsang, C.-F., 2007. Estimating maximum sustainable injection pressure during geological sequestration of CO2 using coupled fluid flow and geomechanical fault-slip analysis. Energy Conv. Man., 47:1798-1807.

Rutqvist J., Rinaldi, A.P., Cappa, F., and Moridis G.J. 2013. Modeling of fault reactivation and induced seismicity during hydraulic fracturing of shale-gas reservoirs. Journal of Petroleum Science and Technology, 107, 31–44 (2013).

Samuelson, J., Spiers, C.J., 2012. Fault friction and slip stability not affected by Co2 storage: Evidence from short-term laboratory experiments on North Sea reservoir sandstones and caprocks. Int. J. Greenhouse Gas Control (in press) dx.doi.org/10.1016/j.ijggc.2012.09.018.

Warpinski N.R. 2009. Microseismic Monitoring: Inside and Out. Journal of Petroleum Technology, November 2009. Society of Petroleum Engineers, SPE 118537.

Warpinski, N.R., M.J. Mayerhofer, K. Agarwal, and J. Du. 2012. Hydraulic fracture geomechanics and microseismic source mechanisms. Society of Petroleum Engineers J. 17:766–780.

Yang, Y., Zoback M.D., 2014. The role of preexisting fractures and faults during multistage hydraulic fracturing in the Bakken Formation. Interpretation. 2, SG25-SG39.

Yong S., Kaiser P.K., Loew S., 2010. Influence of tectonic shears on tunnel-induced fracturing. International Journal of Rock Mechanics & Mining Sciences, 47, 894–907.

Zoback, M.D. 2007. Reservoir Geomechanics. Cambridge University Press, Cambridge, UK.

Zoback, M., Kitasei, S., Copithorne, B., 2010. Addressing the Environmental Risks from Shale Gas Development Briefing Paper1. Natural Gas and Sustainable Energy Initiative, World watch Institute, Washington, July2010 (http://www.worldwatch.org).
<table><caption>Table 1. Rock characteristics considered in the base-case simulation</caption><thead><tr><td>Parameters</td><td>Shale</td><td>Fault</td></tr></thead><tbody><tr><td>Young's modulus, E (GPa)</td><td>30</td><td>30</td></tr><tr><td>Poisson's ratio, v (-)</td><td>0.2</td><td>0.2</td></tr><tr><td>Rock density, ρ<sub>s</sub> (kg/m<sup>3</sup>)</td><td>2700</td><td>2700</td></tr><tr><td>Biot's coefficient, α (-)</td><td>1</td><td>1</td></tr><tr><td>Matrix friction angle, φ (°)</td><td>-</td><td>75</td></tr><tr><td>Matrix cohesion (MPa)</td><td>-</td><td>6</td></tr><tr><td>Matrix tensile strength (MPa)</td><td>-</td><td>1.6</td></tr><tr><td>Joint peak friction angle, φ (°)</td><td>-</td><td>31</td></tr><tr><td>Joint residual friction angle, φ (°)</td><td>-</td><td>22</td></tr><tr><td>Joint cohesion (MPa)</td><td>-</td><td>3</td></tr><tr><td>Joint residual cohesion</td><td>-</td><td>0</td></tr><tr><td>Joint tensile strength (MPa)</td><td>-</td><td>1.6</td></tr><tr><td>Joint dilation angle, ψ (°)</td><td>-</td><td>10</td></tr><tr><td>Permeability, k (m<sup>2</sup>)</td><td>10<sup>-19</sup></td><td>10<sup>-19</sup></td></tr></tbody></table>
Figure captions

1

2

3

4

5

6

7

Figure 1. (Left) 2D-view of the boundary conditions used in the (Right) 3D model geometry for simulating injection and hydraulic fracturing within an 80° dipping fault.

Figure 2. Results of a base-case simulation for a 3-hour long hydraulic fracturing stage considering injection at 2000 m depth into a fault dipping 80 degrees, and having a peak coefficient of friction of 0.6 and a residual friction of 0.4: (a) Downhole well pressure and (b) fault slip and events moment magnitude as function of time, and (c to e) pore pressure, plastic shear strain (signifying shear slip), and plastic tensile strength (signifying tensile failure) at the end of the 3-hour injection.

8

9

10

11

12

13

14

Figure 3. Sensitivity of simulation results to changes in injection depth: Vertical extent of (a) plastic tensile strain (signifying tensile failure) and (b) plastic shear strain (signifying shear slip), and time evolution of (c) well pressure and (d) slip magnitude, and (e) magnitude distribution for each value of injection depth and (f) frequency and maximum magnitude as a function of the injection depth.

15

16

17

18

19

20

Figure 4. Sensitivity of simulation results to changes in fault dip: Vertical extent of (a) plastic tensile strain (signifying tensile failure) and (b) plastic shear strain (signifying shear slip), and time evolution of (c) well pressure and (d) slip magnitude, and (e) magnitude distribution for each dip angle and (f) frequency and maximum magnitude as a function of fault dip angle.

21

22

23

24

25

26

Figure 5. Sensitivity of simulation results to changes in the critical slip distance (or equivalent plastic shear strain threshold) for shear strength drop: Vertical extent of (a) plastic tensile strain (signifying tensile failure) and (b) plastic shear strain (signifying shear slip), and time evolution of (c) well pressure and (d) slip magnitude, and (e) magnitude distribution for each value of the plastic shear strain threshold and (f) frequency and maximum magnitude as a function of plastic shear strain threshold.

27

28

29

30

31

32

33

Figure 6. Sensitivity of simulation results to changes in residual coefficient of friction: Vertical extent of (a) plastic tensile strain (signifying tensile failure) and (b) plastic shear strain (signifying shear slip), and time evolution of (c) well pressure and (d) slip magnitude, and (e) magnitude distribution for each value of the residual friction and (f) frequency and maximum magnitude as a function of residual coefficient of friction.

34

35

36

37

38

39

40
<figure><img src="image_2.png" /><figcaption>Figures
Figure1</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure2</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure4</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure5</figcaption></figure>