Simulation of Injection and Production and MEQ in ... - Semantic Scholar

Simulation of Injection and Production and MEQ in Large-Scale Fracture Networks Moien Farmahini-Farahani, Ahmad Ghassemi Mewbourne School of Petroleum & Geological Engineering, the University of Oklahoma, Norman, OK, 73019, U.S.A.

Abstract Economical production from engineered fractures of reservoirs requires reliable estimation of fluid flow, permeability, and triggered micro-seismicity in large fracture networks. In this work a fast multipole displacement discontinuity method is employed to simulate the response of large-scale natural fracture networks to injection and production. This simulation includes pressure change inside the fracture network, induced normal and shear displacement, slip, dilation, friction coefficient, seismic moment, and micro-earthquakes under anisotropic in-situ stresses acting on a randomly generated set of fractures. Rate-and-state friction model and Mohr-Coulomb criteria are used to model slip and micro-earthquakes. Case studies, representing two fracture network scales having different fracture distributions are presented. Results show that increase in pressure inside fully connected fractures leads to slip events and nucleation of high magnitude micro-earthquakes. Orientation of fractures and fluid pressure interaction are influential factors for nucleation of earthquake slip events. Keywords: Large scale fracture network, Fast multipole method, Displacement discontinuity method, Rate-and-state friction, Micro-earthquakes. 1. Introduction The deformation of natural fracture networks is a crucial aspect of hydraulic stimulation of unconventional resources. Estimating permeability and fluid circulation in a large-scale fracture network in rock mass as a result of injection/production can benefit engineering fractured reservoirs. Slip on fractures increases reservoir permeability through dilatational opening and also, results in micro-seismicity during the injection and production. Changes in the flow pressure and induced deformation in the fracture networks might cause a joint to separate or slip which represents a micro-earthquake (MEQ) occurrence. Realistic models taking into account the shear deformation of mechanically-interacting fractures can adequately represent a fractured rock. These features can be effectively captured using the Displacement Discontinuity Method. The Displacement Discontinuity Method (DDM) has some advantages over other methods, such as the finite element method. For the DDM, instead of meshing the whole field, only fracture surfaces (boundary) are meshed. This means lower number of unknowns or lower number of equations that need to be solved. However, since the influence of all fracture elements is required to be calculated, the equation system matrix is non-symmetric and fully populated. In other words, the conventional DDM is not computational efficient for large fracture systems. Fast approximation techniques, such as the Fast Multipole Method (FMM) (Liu, 2009, Liu and Nishimura, 2006), Hierarchical Matrices (H-Matrices) (Benedetti et al., 2008), wavelet method (Beylkin et al., 1991), panel clustering (Aimi et al., 2005), or fast Fourier transformation methods 1

(Phillips and White, 1997) can significantly improve the computational cost and required memory for boundary element problems. The two most common approximation methods are FMM and HMatrices. Comparison of these two methods shows that FMM requires less memory and builds the coefficient matrix faster and once for the entire simulation (Brunner et al., 2010, Forster et al., 2003, Buchau et al., 2003). Thus, FMM is potentially a more desirable approximation method to be combined with DDM for simulation of large-scale fracture networks. Recently, pressure-displacement coupled FMDDM has been successfully applied to accelerate the calculation for fracture networks in reservoir geomechanics problems (Verde and Ghassemi, 2013a, b, Farmahini-Farahani and Ghassemi, 2015). Presented results in these works are promising because response of the fracture network to injection and production can be accurately captured for regular and natural fracture networks. However, in a coupled FMDDM approach, the coefficient matrix has some degree of inconsistency and is larger compared to a sequential approach. Therefore, the solution requires a higher number of iterations to converge and it is slower compared to the sequential method. In a sequential FMDDM approach, pressure inside the fractures is calculated first and then displacement is calculated based on the change in the pressure. In this work, we use a sequential FMDDM to initially calculate stress and displacement of largescale natural fracture networks and then estimate, based on calculated slip, induced mircoearthquakes on the fracture surfaces. Also, in this study, change in state of elements, as a result of fracture deformation, from stick to slip or separation (opening), is taken into account. Mohrcoulomb criterion along with rate-and-state friction are employed to examine slip events and the friction coefficient. The FMDDM is verified by comparing results with an analytical solution and a conventional DDM code. Simulated cases with 56000 degrees of freedom are presented to highlight the adaptability of the FMDDM method, hydraulic natural fracture deformation, and induced micro-earthquakes. Hopefully the application of approximate method, such as the Fast Multipole Method (FMM) combined with DDM, with capability of simulating large-scale fracture networks in rock masses, will lead to a better understanding of the response of the reservoir to injection and production, as well as more effective techniques to design reservoirs. 2. Numerical Model 2.1. Sequential pressure and displacement Injection or/and production lead to change of pressure in fluid flow inside a fracture network and consequently, change of stress over the fracture surface and aperture-size change follow. To calculate flow through the interconnected fractures, mass balance and Darcy’s law are combined as: 𝜕𝑞 𝜕𝑝 (1) = 𝑤𝑓 𝐿𝑐𝑡 − 𝑞𝑖 𝜕𝑠 𝜕𝑡 where, 𝑤𝑓 is fracture aperture, 𝐿 is element length, 𝑐𝑡 is total compressibility, and 𝑞𝑖 is production 𝜕𝑞

rate per thickness. In Eq. (1), 𝜕𝑠 is net flow rate inside a fracture, and Darcy’s law indicates that 𝑞 = −

𝑤𝑓 𝑘 𝜕𝑝 𝜇

𝜕𝑠

, where, 𝑘 is fracture permeability and 𝜇 is fluid viscosity. The first term on the right

represents change of fluid volume as a result of fracture deformation and expansion/compression. 𝑤𝑓 𝑐𝑡 is inverse of the total stiffness of the fracture surface and the fluid which are acting like a series of two springs with stiffness of 𝐾𝑎 = 𝐾𝑛 ⁄(1 − 𝑛) and 𝐾𝑓 = 1/𝑤𝑓 𝑐𝑡 , respectively, where, 2

𝑐𝑜𝑛𝑡𝑎𝑐𝑡 𝑎𝑟𝑒𝑎

𝐾𝑛 is normal stiffness, 𝑛 is contact porosity 𝑛 = 1 − 𝑗𝑜𝑖𝑛𝑡 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 𝑎𝑟𝑒𝑎 , and 𝑐𝑓 is fluid compressibility. 𝑠 represents fluid flow direction which depends on the connectivity of fractures.

First, the flow pressure is calculated from eq. (1) and then the induced displacement, as a result of change in pressure, is calculated. Induced displacement is solved through Fast Multipole Displacement Discontinuity Method (FMDDM), in which the Fast Multipole Method (FMM) is incorporated with Displacement Discontinuity Method (DDM) to accelerate the solution for largescale problems. The DDM is an indirect boundary element method suitable for analysis of fracture deformation (change in aperture size) in a medium containing multiple fractures (Crouch and Starfield, 1983). It is based on the analytical solution to the problem of a finite line segment centered in the x-y plane with constant normal and shear displacement discontinuities (See figure 1a).

Figure 1. (a) A two-dimensional fracture segment with constant discontinuity displacement components D x and Dy; (b) Elemental segments over a fracture.

For calculation of the displacement over a fracture, the fracture is divided into N elemental segments (see figure 1b). Induced normal and shear stresses on the i-th fracture segment is calculated as the sum of effects of all N elements (Crouch and Starfield, 1983). For natural fractures with in-situ stress, the fluid flow pressure and fracture stiffness set of equations can be rearranged as: 𝑁

0=

𝑁

𝑖𝑗 𝑗 ∑ 𝐴𝑠𝑠 𝐷𝑠 𝑗=1

𝑖𝑗

𝑗

+ ∑ 𝐴𝑠𝑛 𝐷𝑛 + 𝐾𝑠𝑖 𝐷𝑠𝑖 𝑗=1

𝑁

(2)

𝑁 𝑖𝑗

𝑗

𝑖𝑗

𝑗

−∆𝑃𝑖 = ∑ 𝐴𝑛𝑠 𝐷𝑠 + ∑ 𝐴𝑛𝑛 𝐷𝑛 + 𝐾𝑛𝑖 (𝐷𝑛𝑖 + tan 𝜙𝑑 𝐷𝑠𝑖 ) 𝑗=1 𝑖𝑗

𝑗=1 𝑖𝑗

𝑖𝑗

𝑖𝑗

where, 𝐴𝑛𝑠 , 𝐴𝑛𝑛 , 𝐴𝑠𝑠 , and 𝐴𝑛𝑠 are boundary influence coefficients and 𝐾𝑠 and 𝐾𝑛 are the crack shear and normal stiffness, respectively. Normal stiffness is estimated by a hyperbolic equation as a function of the maximum closure and initial stiffness (Bandis et al., 1983) and 𝜙𝑑 is dilation angle which accounts for change of the normal displacement by change in shear displacement. In this paper, we use the approach employed by (Zhou and Ghassemi, 2011) to incorporate the flow pressure, displacement, stress, and finally estimate the MEQ. Change in fracture permeability because of shear displacement is taken into account. Aperture change for each element is expressed as: 3

∆𝐷𝑛 = 𝐷𝑛 + 𝐷𝑑𝑖𝑙 = (

−∆𝜎𝑛′ ) + 𝐷𝑠 tan(𝜙𝑑 ) 𝐾𝑛

(3)

where, ∆𝜎𝑛′ is the increment of the normal effective stress, 𝐷𝑑𝑖𝑙 is the dilation induced aperture change generated by slip. For pressurized fracture modeling, each element can be in three different states; “slip”, “stick”, or “separation”. Since state of each element might change from one time step to another step, its state must be determined at every time step. An element can be open or closed under the applied fluid pressure and in-situ stress. To determine the state of each element, first the acting load needs to be calculated. Under different state conditions, different left-hand sides (known parameters) for eq. (2) should be considered to calculate displacement discontinuity. If the difference of pressure and normal stress acting on the surface of an element is less than the product of cohesion and cotangent of effective friction angle, the element is open under the load and its state is “separation”: 𝜎𝑛′ = 𝑐 × cot(𝜙𝑒𝑓𝑓 ) = 𝑐 × cot(𝜙𝑓 + 𝜙𝑑 )

(4)

where, 𝜎𝑛′ is the effective normal stress, 𝑐 is cohesion, and 𝜙𝑒𝑓𝑓 is the effective friction angle element surface which is the sum of friction angle and dilation angle. The next step is to determine whether a closed element is in a state of “stick” or “slip”. For this purpose, the Mohr-Coulomb (MC) criterion (Eq. (5)) is used: |𝜎𝑠 | > ±𝜏 = 𝑐 + 𝜎𝑛′ tan(𝜙𝑒𝑓𝑓 )

(5)

where, 𝜏 is the yield stress. If the shear stress of an element exceeds the Mohr-Coulomb criterion, the element is yielding and its state changes to “slip”. To calculate the displacement for such an element, its 𝐾𝑠 is set to zero and the acting shear stress on the element is set to ±𝜏 − 𝜎𝑠0 , where, 𝜎𝑠0 is the in-situ shear stress acting on the surface of the element. At this stage, since the acting shear stress has changed, displacements and then stresses should be recalculated, followed by the Mohr-Coulomb criterion check. This process is long and might need to be repeated several times to reach consistency in state, displacements, and stresses. Normal displacement might slightly change, so we can assume that the normal displacement does not change and the only recalculate the shear displacement through a correlation (Jing et al., 2000): 𝐷𝑠 =

𝑅(𝜏 − 𝜎𝑛′ tan(𝜙𝑒𝑓𝑓 )) 𝛾𝐺

(6)

where, 𝛾 is a geometric parameter and its value depends on shape of the fracture (Dieterich, 1992), 𝐺 is shear modulus of the encompassing rock, and 𝑅 is the fracture half length. When small rates of displacement are considered, Eq. (6) is in good agreement with (Barton et al., 1985). When an element is not surpassing the Mohr-Coulomb criterion, its state is “stick”. In this state, normal stiffness is estimated through the equation proposed by Barton-Bandis, Eq. (7) (Bandis et al., 1983) and the shear stiffness is a constant equaling the initial value. −2

𝜎𝑛′ ] 𝐾𝑛 = 𝐾𝑛𝑖 [1 − 𝑉𝑚 𝐾𝑛𝑖 + 𝜎𝑛′

(7)

where, 𝐾𝑛𝑖 is initial normal stiffness and 𝑉𝑚 is the maximum closure. 4

2.2. Fault friction and micro-earthquakes Rate-and-state friction predicts friction coefficients based on the slip rate 𝑣 and the single-state variable 𝜃 which represents the slip history of fracture. The friction coefficient according to the rate-and-state friction is (Segall, 2010): 𝑓 = 𝑓0 + 𝑎𝑙𝑛

𝑣 𝜃 + 𝑏𝑙𝑛 𝑣0 𝜃0

(8)

where, 𝑎 and 𝑏 are constants measured experimentally, 𝑣 is the slip velocity, 𝑓0 is the initial coefficient such that 𝑓 = 𝑓0 when 𝑣 = 𝑣0 and 𝜃 = 𝜃0 , and 𝜃 is a single-state variable calculated by aging law (Eq. (9) below) reflecting the sliding history effects. Since contacts are destroyed and built under slip, the state variable presents a measure of loading support provided by the sliding surface and contacts strengthen with age (Ruina, 1980): 𝑑𝜃 𝑣𝜃 = 1− 𝑑𝑡 𝑑𝑐

(9)

𝑑𝜃

When 𝑑𝑡 = 0, the characteristic slip distance is 𝑑𝑐 = 𝜃0 𝑣0 . The magnitude of characteristic slip distance depends on the surface roughness and fault gouge (Dieterich, 1979). Here a constant slip velocity is assumed. The seismic moment measures the size of energy release of an earthquake caused by a slip displacement and it is estimated as an integral of the shear modulus times the shear displacement over the fracture area that slip: 𝑀0 =

∫

𝐺𝑆 𝑑𝐴

(10)

𝑆𝑙𝑖𝑝 𝑎𝑟𝑒𝑎

where, G is the shear modulus of the rock, 𝑆 is cumulative slip, and 𝐴 is the slip area. The magnitude of earthquake produced by the seismic moment can be estimated as (McGarr et al., 1979): 2 𝑀 = log 𝑀0 − 10.7 3 where 𝑀0 is in dyn/cm.

(11)

3. Results and Discussion The numerical simulation presented here makes it possible to analyze the response of fracture networks to fluid injection and production. In this study fluid is injected and produced at the same rate from two wells (one for the injection and one for the production). Change in the normal and shear displacement, conductivity, pressure, friction angle and induced micro-seismicity have been examined to quantify affected fractures. These parameters indicate how the fracture permeability is affected in the rock mass as a result of injection and production. In this study, fracture propagation is ignored.

5

3.1. Sample case 1- A fully connected fracture network Figure 2 shows the geometry of a fracture network consisting of 100 fully connected fractures. The model represents a two-dimensional horizontal section of a reservoir with a unit thickness. It is assumed that injection and production wells are vertical wells located as marked in figure 2. The fracture network expands over a rock mass of size 1100 m × 700 m. The state of stress in the plane of the model is anisotropic, with a higher in the Y-direction. Applied injection and production rates are per thickness of the model. It is assumed that the fractures all have a uniform aperture at the beginning and are conductive.

Figure 2. Geometry of the fracture network. Injection and production wells are marked with red circles.

Initial variables and geomechanical parameters are listed in Table 1. Table 1. Assumed parameters and initial values of variables.

Parameter In-situ 𝜎𝑥𝑥

Value 11 MPa

Parameter Viscosity

In-situ 𝜎𝑦𝑦 Initial pressure Poisson’s ratio Shear module Initial aperture size Shear stiffness Normal stiffness

13 MPa

Injection rate (per thickness)

0.001 m /s

Production rate (per thickness)

0.001 m /s 2° 30°

Dilation angle Intact friction angle

Value 2.03×10

Compressibility

-10 -3

MPa.s -1

0.42×10 MPa Simulation time 1000 s Time steps 10 Rate-and-state friction parameters a 0.0110 b 0.0125 -6 dc 6×10 m

10 MPa 0.25 15 GPa 1 mm 50 GPa/m 25 GPa/m 3

V0

3

-9

5×10 m/s

This is a static model without fracture propagation and constant in-situ stresses. The fracture aperture size is scaled up to 3000 times larger than the actual size to be easily displayed; otherwise, millimeter-sized fractures cannot be seen over such large areas. Figure 3a shows change in pressure through the fracture network as a result of the injection and 6

production. With continued operation, pressure reaches equilibrium through the fracture network. The pressure around the injection well is higher and leads to lower effective normal stress around the higher-pressure zone. Pressure is lower around the production well, which causes higher effective normal stress on fracture surfaces in the lower-pressure zone. Since, the state of stress is anisotropic, different fracture orientations respond differently because of different initial normal and shear stresses acting on the fractures, resulting in potentially different levels of slip, and dilatation (see Figure 3b).

Figure 3. (a) Pressure distribution inside the fracture network; (b) Effective normal stress.

Figure 3b clearly shows that fractures in the higher pressure zone experience lower effective pressure. Because the in-situ stress component 𝜎𝑦𝑦 is higher than the 𝜎𝑥𝑥 component, the fractures with higher inclination from the x-axis tend to have a lower effective normal stress. The same situation applies to the low-pressure zone, i.e., a higher inclination leads to a lower effective normal stress. The lower effective normal stress reduces the shear strength and increases the possibility of failure according to the Mohr-Coulomb criterion (Eq. (5)). Note that since the pressure inside one fracture is nearly uniform, when one element of a fracture meets the MC criterion, it is highly possible that all other elements of the fracture would also fail; thus a fracture is said to fail when all of its elements fail. In Figure 4a, fractures that met the MC are shown in red. Note that all fractures in red are located in the high-pressure zone and the effective normal stress acting on their surfaces is less than 2 MPa. When an element fails, it slides with an associated dilatational opening (for slip less than 3mm). Figures 4b and 4c show fracture slip and related dilatational opening. Elements with higher slip have higher dilation (see Figure 4c). Figure 4c shows dilation of fracture surfaces as a result of shear displacement. Fracture dilation associated with slip increases fracture aperture and permeability.

7

Figure 4. (a) Status of fractures based on the MC criterion; (b) Slip over fracture surfaces; (c) Dilatational opening.

In addition, slip on the fracture surfaces potentially induces micro-earthquakes (Figure 5). For each fracture element, the hypocenter is defined at its center. For a fault/fracture, the MEQ is obtained through summation of individual fracture element slip values. The estimated seismic moment is per unit fracture thickness (depth). Higher magnitude MEQs occur around the higher pressure zone where the effective normal stress is lower, and slipping and dilatational opening are larger.

Figure 5. Induced MEQs hypocenter map. Sized by seismic moment.

3.2. Sample case 2- A partially connected fracture network The fast multipole displacement discontinuity method is computationally efficient compared to the conventional DDM. This is important especially when the fracture network is a large-scale and

8

complex natural fracture network, as is the case of Enhanced Geothermal System (EGS). A complex large-scale fracture network consisting of 700 fractures is presented in Figure 6.

Figure 6. Geometry of the fracture network; inset shows zoom-in view for injection and production wells.

Figure 6 shows a 2D horizontal section of the large-scale fracture network in a 5000 m ×5000 m domain, with 200 fully connected fractures and 500 partially connected or single fractures around the main network. The main network comprises the fully connected fractures and is surrounded by other fractures (figure 6 inset). Each fracture is divided into several elements and the total number of elements is 28,334 (56,668 degrees of freedoms). Injection and production wells are perpendicular to horizontal section and are marked in the inset. The state of stress is anisotropic, with a higher in-situ stress component in the y-direction. It is assumed that all fractures are conductive and have uniform apertures at the beginning of simulation. All initial variables and geomechanical parameters are the same as in sample case 1 (Table 1). Figure 7a and 7b show the pressure distribution inside the fracture network. Pressure is higher and lower than the initial pore pressure around the injection well and the production well, respectively. Figure 7b provides a zoom-in view, in the main fracture network, around the injection and production wells. Injection and production rates are not high for this geometry and pressure does not change significantly inside the fracture network. However, even small change in pressure (~150 KPa) induces normal displacement and slip events over the fracture surfaces. The role of fractures orientations is evident. Fractures with higher inclination from the x-axis are prone to fail the MC criterion as a result of small increase in the pressure. Figure 7c shows change in aperture as a result of the injection and production. Higher pressure increases the fracture apertures around the injection well and, likewise, fracture apertures decrease around the production well.

9

Figure 7. (a) Pressure distribution inside the large-scale fracture network; (b) Zoom-in view of the pressure distribution around the injection and production wells; (c) Change in the aperture around the injection and production wells.

Figure 8a shows slip on the fracture surfaces in the main fracture network. An increase in the pressure causes a decrease in effective normal stress and leads to slip on the fracture surfaces, subsequently induces MEQs. Induced MEQs hypocenter map as a result of injection and production for the fracture network and corresponding slip are shown in figure 8b. MEQ occurrence is sized by seismic moment measured over the slip area. Three high magnitude MEQs occur in the main network correspond to slip occurrences shown in figure 8a. These events are directly induced by the decrease in effective normal stress. This leads to the MC failure and occurrence of slip over fracture surfaces, which are orientated in the direction of the shear stress acting on their surfaces. Induced field stress 𝜎𝑦𝑦 (compression positive) is shown in Figure 9. We can see that tensile stresses are created at the fracture tips and compressive stresses are created along the fractures and between fractures. Zoom-in views clearly show that at the fracture tip kidney-like stress distributions are formed.

10

Figure 8. (a) Slip on the main fracture network; (b) Induced MEQs hypocenter map.

Figure 9. Induced field stress distribution in Y direction.

11

4. Conclusions The fast multipole displacement discontinuity method (FMDDM) is employed to analyze the response of two-dimensional fracture networks to injection and production. The sample cases presented in this study are representative of large-scale fracture networks with anisotropic in-situ stress and randomly orientated fully connected fractures. In-situ conditions, such as the magnitude of field stresses and fracture orientations play essential roles in shear stress acting on fracture surfaces. Increase in pressure inside fractures reduces effective normal stress and favor slip deformation subsequently enhancing permeability through dilatational opening and induced MEQs. Higher magnitude MEQs occur around the high-pressure zone in the main network in which fluid is injected and produced. The results presented here show that FMDDM can be used to model large-scale and complex fracture networks. Also, a comprehensive approach, which takes into account the interaction of all elements, shear strength failure, and induced MEQs, is required for a more reliable reservoir design and to better predicts the response of large-scale fracture networks. References Aimi, A., M. Diligenti, and F. Lunardini, 2005. "Panel clustering method and restriction matrices for symmetric Galerkin BEM." Numerical Algorithms, v. 40 (4), p. 355-382. Bandis, S.C., A.C. Lumsden, and N.R. Barton, 1983. "Fundamentals of rock joint deformation." International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, v. 20 (6), p. 249-268. Barton, N., S. Bandis, and K. Bakhtar, 1985. "Strength, deformation and conductivity coupling of rock joints." International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, v. 22 (3), p. 121-140. Benedetti, I., M.H. Aliabadi, and G. Davı, 2008. "A fast 3D dual boundary element method based on hierarchical matrices." International Journal of Solids and Structures, v. 45, p. 2355. Beylkin, G., R. Coifman, and V. Rokhlin, 1991. "Fast Wavelet Transforms and Numerical Algorithms." Communications on Pure and Applied Mathematics, v. 44 (2), p. 141-183. Brunner, D., M. Junge, P. Rapp, M. Bebendorf, and L. Gaul, 2010. "Comparison of the fast multipole method with hierarchical matrices for the helmholtz-bem." Computer Modeling in Engineering & Sciences, v. 58 (2), p. 131-158. Buchau, A., W.M. Rucker, O. Rain, V. Rischmüller, S. Kurz, and S. Rjasanow, 2003. "Comparison between different approaches for fast and efficient 3D computations." IEEE Transactions on Magnetics, v. 39 (3), p. 1107-1110. Crouch, S.L., and A.M. Starfield, 1983. "Boundary element methods in solid mechanics: With applications in rock mechanics and geological engineering", 322. London, George Allen & Unwin. Dieterich, J.H., 1979. "Modeling of rock friction: 1. Experimental results and constitutive equations." Journal of Geophysical Research: Solid Earth v. 84 (B5), p. 2161–2168. 12

Dieterich, J.H., 1992. "Earthquake nucleation on faults with rate- and state-dependent strength." Tectonophysics, v. 211, p. 115-134. Farmahini-Farahani, M., and A. Ghassemi, 2015. "Analysis of fracture network response to fluid injection ". Stanford University, California, USA. Forster, H., T. Schrefl, R. Dittrich, W. Scholz, and J. Fidler, 2003. "Fast boundary methous for magnetostatic interactions in micromagnetics." IEEE Transactions on Magnetics, v. 39 (5), p. 2513-2515. Jing, Z., J. Willis-Richards, K. Watanabe, and T. Hashida, 2000. "A three-dimensional stochastic rock mechanics model of engineered geothermal systems in fractured crystalline rock." Journal of Geophysical Research: Solid Earth, v. 105 (B10), p. 23663-23676. Liu, Y., 2009. "Fast Multipole Boundary Element Method", Cambridge University Press. Liu, Y.J., and N. Nishimura, 2006. " The fast multipole boundary element method for potential problems: A tutorial." Engineering Analysis with Boundary Elements, v. 30 (5), p. 371-381. Mcgarr, A., S.M. Spottiswoode, N.C. Gay, and W.D. Ortlepp, 1979. "Observations relevant to seismic driving stress, stress drop, and efficiency." Journal of Geophysical Research: Solid Earth v. 84 (B5), p. 2251-2261. Phillips, J.R., and J.K. White, 1997. "A Precorrected-FFT Method for Electrostatic Analysis of Complicated 3-D Structures." IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, v. 16 (10), p. 1059-1072. Ruina, A. 1980. Friction laws and instabilities: A quasistatic analysis of some dry frictional behavior. Phd Brown University, Providence, RI. Segall, P., 2010. "Earthquake and Volcano Deformation". Princeton, New Jersey, Princeton University Press. Verde, A., and A. Ghassemi, 2013a. "Efficient Solution of Large-Scale Displacement Discontinuity Problems using the Fast Multipole Method." Proc., 47th US Rock Mechanics / Geomechanics Symposium., San Francisco, CA, USA,. Verde, A., and A. Ghassemi, 2013b. "Fracture Network Response to Injection Using An Efficient Displacement Discontinuity Method." Proc., 37th Annual meeting Geothermal Resources Council (GRC), Las Vegas, NV, USA,. Zhou, X., and A. Ghassemi, 2011. "Three-dimensional poroelastic analysis of a pressurized natural fracture." International Journal of Rock Mechanics & Mining Sciences, v. 48, p. 527-534.

13