We investigate fracture propagation induced by hydraulic fracturing with water-injection, using numerical simulation. For full 3D rigorous modeling, we employ a numerical method that can model failure due to tensile and shear stresses, dynamic nonlinear permeability, the dual continuum approach, leak-off in all directions, and thermo-poro-mechanical effects. From the numerical results, fracture propagation is not the same as propagation of the water front, because fracturing is governed by geomechanics, whereas water saturation is determined by fluid flow. At early times, the water saturation front is almost identical to the fracture tip, showing that the fracture is mostly filled with injected water. However, at late times, advance of the water front is retarded, compared to the fracture propagation, yielding a significant gap between the water front and the fracture top, which is filled with reservoir gas. We also find considerable leak-off of water to the reservoir. The inconsistency between the fracture volume and the volume of injected water cannot properly estimate the fracture length, when it is estimated based on the simple assumption that the fracture is fully saturated with injected water. As an example of flow-geomechanical responses, we identify pressure fluctuation under constant water injection, because hydraulic fracturing is itself a set of many failure processes, in which pressure drops every time when failure occurs. The fluctuation decreases as the fracture length grows.

We also study application of electromagnetic (EM) geophysical methods, because the EM geophysical methods are highly sensitive to changes in porosity and pore-fluid properties, such as water injection into gas reservoirs. We employ a 3D finite-element EM geophysical simulator and evaluate the sensitivity of the crosswell EM method for monitoring fluid movements in shaly reservoirs. For the sensitivity evaluation, reservoir models are generated through the coupled flow-geomechanical simulator and are transformed via a rock physics model into electrical conductivity models. It is shown that anomalous conductivity distribution in the resulting models is closely related with injected water saturation but little with newly-created unsaturated fractures. The numerical modeling experiments demonstrate that the crosswell EM method can be highly sensitive to conductivity changes that directly indicate the migration pathways of the injected fluid. Accordingly, the EM method can serve as an effective monitoring tool for distribution of injected water (i.e. migration pathways) during hydraulic fracturing operations.

Hydraulic fracturing with multiple horizontal wells has widely been used in stimulating abundant shale gas reservoirs in order to increase reservoir permeability and productivity (Zoback 2007; Fjear et. al. 2008; Fisher and Warpinski 2012; Zoback et. al. 2010; Cipolla et. al. 2010). Creation of fractures increases permeability of geological formations by several orders, which can meet economical feasibility of gas production of very low permeable reservoirs such as tight gas and shale gas reservoirs. After the success of development in the Barnett shale, other shale gas plays such as the Marcellus, Eagle Ford, Haynesville, New Albany, Antrim, and Fayetteville shales in the United States have been considered as potential resources near in the future. Massive hydraulic fracturing with multiple stage stimulation can significantly enhance reservoir permeability followed by productivity, increasing stimulated rock volume (SRV) (Kargbo et. al. 2010; Arthur et al. 2008). On the other hand, the issues on the environmental impact of hydraulic fracturing have been raised by environmental scientists, such as unstable growth of the fractures, groundwater contamination, and induced seismicity from fault activation. For example, Osborn et al. (2011) claimed from the isotope analysis that methane dissolved in groundwater originated from shale gas reservoirs. Geomechanical failure along the well casing due to incomplete cementing might also be problematic, yielding a potential pathway of methane from the reservoirs to the drinking water aquifers (Zoback et al. 2010). Reservoir engineers however argue that the environmental impact is still limited. For example, according to Fisher and Warpinski (2012), the fracture lengths based on measured microseismic data are finite and the fracture tips are far below the drinking water zones. Rutqvist et al. (2013) investigated induced seismicity by numerical simulation, showing that magnitudes of the induced seismicity by fault activation are far below the level discernible by human.

Much research in simulating fracture propagation has been performed on hydraulic fracturing for several decades in order to estimate SRV and the lengths of the created fractures, mainly based on 2D or pseudo 3D (Perkins and Kern, 1961; Nordren, 1972; Adachi et.al., 2007). Even though those models can provide numerical efficiency and reduce computational cost, full coupling in 3D between flow and geomechanics is required for rigorous modeling of fracture propagation and more reliable risk assessment. Recently, Ji et al. (2009) and Dean and Schmidt (2009) performed numerical modeling of full 3D hydraulic fracturing in the context of coupled flow and geomechanics. Furthermore, Kim and Moridis (2013) extended their modeling, accounting for the dual continuum model as well as simultaneous tensile and shear failure.

Hydraulic fracturing is typically based on water injection, for which water is mixed with some proppants and additional chemicals (King, 2012). Water can easily be obtained and pressurized. However, after hydraulic fracturing operations, it is difficult to withdraw the injected water, which might cause significant decrease of productivity or adverse chemical reactions (Page and Miskimins 2009; Ribeiro and Sarma 2013). Thus, it is critical to accurately estimate fracture propagation as well as fluid flow and proppant migration.

In this paper, we numerically investigate fracture propagation by water injection for hydraulic fracturing operations. We focus on fracture propagation and water movement. The fracture propagation might not be the same as propagation of the water front, because fracturing is governed by geomechanics whereas water saturation is determined by fluid flow. The inconsistency between the fracture volume and the volume of the injected water cannot properly estimate dimension of the fracture when we simply assume the fracture to be fully saturated with water. For full 3D rigorous modeling, we employ a numerical method shown in Kim and Moridis (2013), which can model tensile failure due to normal and shear stresses, dynamic nonlinear permeability and geomechanical moduli, leak-off in all directions during hydraulic fracturing, and thermo-poro-mechanical effects. In this study, we will find a significant gap between the water front and the fracture top, as well as considerable leak-off of water to the reservoir.

Along with 3D rigorous modeling of fracture propagation and fluid flow, it is also important to develop geophysical methods for monitoring migration pathways of injected fluids. Microseismic methods have been widely used in the past decade for estimation of fracture propagation and geometry (Warpinski et al., 2005; Vermylen and Zoback, 2011). In some cases, however, the magnitude of microearthquakes would be too small to be reliably recorded in practice. In addition, the microseismic inversion for determining fracture length and width highly depends on an initial velocity model and can result in significant ambiguity in hydraulic fracture characterization (Johnston and Shrallow, 2011). To reduce the ambiguity in the microseismic monitoring, one can consider electromagnetic (EM) geophysical methods. When hydraulic fracturing operations force a mixture of fluids and proppant into the formation under high pressure, they generate highly-localized changes in porosity and pore fluids. Because the EM methods are sensitive to porosity, fluid saturation and chemistry of the fluid in the pore space, they are a promising tool for illuminating migration pathways of the injected fluids and proppant and can complement microseismic methods. In this paper, we investigate the sensitivity of an EM geophysical method to electrical conductivity changes that are directly correlated with hydraulic fracturing operations.

Accordingly, this paper consists of two major parts: coupled flow-geomechanical modeling and EM geophysical modeling. The reminder of this paper is organized as follows. First, we briefly describe a T+M coupled flow-geomechanical modeling algorithm and simulate fracture propagation and coupled flow-thermal-geomechanical processes in a tight shale gas reservoir. From the simulation, we also obtain changes in reservoir porosity and permeability, pressure and saturation of gas and water, and the fracture opening and displacement over the domain. The resulting reservoir parameters are employed to construct a series of realistic 3D electrical conductivity models through a rock physics model. Subsequently, the sensitivity of an EM geophysical method to conductivity changes due to hydraulic fracturing operations is evaluated. It is demonstrated that by using an electrically-engineered high-conductivity fluid, we can significantly improve the sensitivity of the EM method to migration pathways of the injection fluid.

Water injection

Two phase flow within a fracture occurs during water injection of hydraulic fracturing processes. As water flows through the fracture, reservoir gas is pressurized while some of water leaks off into the reservoir formation, as shown Fig. 1. In the left figure, a stimulated zone consists of a main fracture, several fissures, and partially continuous small fractures, depending on complexity of hydraulic fractures (Fisher and Warpinski, 2012). We can model flow and geomechanics of different fractures, using the dual continuum (also called double porosity) approach.

Fig. 1

Schematics of complex fractures and fluid movement around the fractures. Hydraulic fractures create different types of the stimulated zone shown in the left figure. In the right figure, multiphase flow depends on viscous force, buoyancy, and poromechanical effects, opening and closing the fractures.

Fig. 1

Schematics of complex fractures and fluid movement around the fractures. Hydraulic fractures create different types of the stimulated zone shown in the left figure. In the right figure, multiphase flow depends on viscous force, buoyancy, and poromechanical effects, opening and closing the fractures.

Close modal

Let us consider the time when fracturing occurs, while there is little pressure diffusion. At this time scale, after the fracturing, pore pressure is still below pore-pressure that can open the created fracture. We can then assume that reservoir gas fills fractured space immediately once fracturing occurs, based on high mobility of gas (i.e., low viscosity and high relative permeability), shown in the right of Fig. 1. When a fracture is closed, water cannot move to the new created fracture. After the fracture becomes open, reservoir gas and the injected water come into the new fracture space by viscous force as well as buoyancy. Furthermore, the increase of the fracture volume can decrease pore-pressure because of poro-mechanical effect, which might allow the reservoir gas to leak into the fracture, if the pore-pressure is below gas pressure of the reservoir.

Consider the case that reservoir pressure and temperature are 17.1MPa and 58.8°C, respectively, which come from a generalized model of Marcellus shale. Viscosities of liquid water and gas phases are 4.79× 10−4Pa·s and 1.78× 10−5Pa·s, respectively. Density of gas at the reservoir condition, 1.15 ×102kg·m−3, cannot also be neglected, compared with density of water, 9.89 × 102 kg·m−3. From the densities and viscosities of water and gas, we can estimate the pattern of water saturation movement of buoyancy, constructing the following equation of the fractional flow curve, fw, after ignoring the gradient of capillary pressure within a fracture.

fw=1Ngkr,g1+kr,gμgkr,wμw,kr,J=max{0,min{(SJSr,J1.0Sr,w)nk,1}},
(1)

where Ng, kr,J, and μJ are the gravity number, relative permeability and viscosity of phase J. The subscript g and w indicate gas and water phases, respectively. SJ and Sr,J are saturation and residual (irreducible) saturation of phase J. We use the modified versions of Stone's relative permeability model.

Fig. 2 shows the fractional flow curves with various gravity numbers. From Fig. 2, the movement of water saturation is estimated to be piston-like displacement, generating a shock wave. On the other hand, an aperture distribution along the fracture from geomechanics is continuous. Reservoir gas would then be accumulated on the upper area of the fracture, which yields a dry zone. We can anticipate that the gas volume within the fracture become larger as the simulated zone becomes large, because reservoir gas that fills the simulated zone moves upward due to buoyancy. There is an additional issue how much water infiltrates into the reservoir rock matrix, because leak-off of water into rock matrix can significantly affect mobility of gas, reducing relative permeability, when gas is produced. Thus, rigorous modeling that can simulate coupled flow and geomechanics is required to predict physical processes during hydraulic fracturing operations accurately.

Fig. 2

Left: fractional flow curves at different gravity numbers. Right: Schematics of the fracture opening and water saturation distribution.

Fig. 2

Left: fractional flow curves at different gravity numbers. Right: Schematics of the fracture opening and water saturation distribution.

Close modal

Failure condition and permeability

Hydraulic fracturing is based on tensile failure, which creates a fracture and opens it. We employ a stress-based failure condition for large scale fracture propagation, rather than a toughness-based condition, and use a failure condition (Ruiz et al., 2000), written as

σc(=β2(tt2+ts2)+tn2)Tc,
(2)

where σ′c is the effective stress for tensile failure, which consists of shear and normal effective stresses. tn, tt and ts are normal and shear effective stresses acting on a fracture plane (Fig. 3), and Tc is tensile strength. Using β of Eq. 2, we account for contribution of shear effective stress to tensile failure.

Fig. 3

Schematics of tensile failure in 3D. Left: planar fracture propagation along the vertical direction. Right: modeling of vertical fracture propagation using no horizontal displacement condition at the plane that contains the vertical fracture (Kim and Moridis, 2013).

Fig. 3

Schematics of tensile failure in 3D. Left: planar fracture propagation along the vertical direction. Right: modeling of vertical fracture propagation using no horizontal displacement condition at the plane that contains the vertical fracture (Kim and Moridis, 2013).

Close modal

When implementing the tensile failure, we implicitly employ a node splitting scheme, as shown in the left of Fig. 3. Then, considering vertical fracture propagation, we can simplify the node splitting scheme into a method that updates outer boundary condition (i.e., from the Dirichlet condition to the Neumann condition, as shown in the right of Fig.3), not introducing an internal boundary. This scheme was previously used in Ji et al. (2009). We then can reduce computational resources and code management effort.

For the modeling of the fracture permeability corresponding to tensile failure, we use the modified Cubic law (e.g., Snow, 1965; Rutqvist and Stephansson, 2003), written as

Qw=acωfnp12μwH(pρwg),
(3)

where ωf is the aperture, also called fracture opening. Qw is flow rate of water. H is the fracture plate width, and g is the gravity vector. For np = 3.0, Eq. 3 reduces to the original cubic law. ac is the correction factor the reflects the fracture roughness. We use a minimum fracture permeability, still much higher than the permeability of intact rock, even when the fracture is closed.

Coupled flow and geomechanics

We employ quasi-static mechanics, based on momentum balance, when solving geological deformation and failure. Governing equations of fluid and heat flow come from conservation laws of fluid mass and energy. In hydraulic fracturing operations, flow and geomechanics are tightly coupled because permeability is a strong function of material failure and deformation of a fracture, as described in Eq. 3. High deformability of fractures also affects coupling in pore volume, which cannot be neglected. The coupling of pore volume can be modeled by poromechanics (Coussy, 2004). Specifically, according to Berryman (2002) and Kim et al. (2012), the porosity coupling for the multiple continuum approach that can represent a fracture-rock matrix system (or a system of the main fracture and small local fractures) can be modeled as

Φl=(αl2Kl+αlΦlKs)δpl+3αT,lαlδTlblηlδσv,bl=αlηlKl,
(4)

where Φ is the Lagrange porosity, defined as the ratio of the pore volume in the deformed configuration to the bulk volume in the reference (initial) configuration. The subscription l indicates subelements in a gridblock. T is temperature. Kl, Ks,ηl, αT,l αl, are the drained bulk modulus, intrinsic solid grain modulus, the volume fraction, the thermal dilation coefficient, and Biot's coefficient for subelement l within a gridblock, respectively. σv is total mean (volumetric) stress at the gridblock.

In numerical modeling, we employ the finite volume method for flow, in which flow variables take piecewise-constant interpolation. For geomechanics, we use the finite element method, taking linear interpolation of displacement. In time discretization, we employ the backward Euler method. These discretizations for flow and geomechanics are widely used in reservoir simulation and computational geomechanics communities, respectively.

We take a sequential method to solve coupled flow and geomechanics problems. Specifically, we use the fixed-stress sequential method in solving two-way coupling in pore-volume, written as

Φln+1Φln=(αl2Kl+αlΦlnKs)Φlncp(pln+1pln)+3αT,lαlΦlncT(Tln+1Tln)blηl(σvnσvn1)ΔΦln,
(5)

where the superscript n is the time level in time discretization. c and cT correspond to the pore-compressibility and thermal expansivity used in reservoir simulation, respectively. ΔΦln is called porosity correction, calculated from the previous geomechanics solutions, in order to capture poromechanical effects.

Permeability is calculated based on explicit treatment, depending on failure status and aperture of the fracture at the previous time step. Permeability is a strong function of geomechanical failure, and thus we take a small time strep size that ensures no fracturing between two fracturing events for further accuracy.

From the given numerical schemes, we can make use of existing flow and geomechanics codes, by constructing an interface between them. In this study, we couple TOUGH+RealGasH2O, a flow simulator, to ROCMECH, a geomechanics simulator, shortly called T+M. T+M takes the several verification tests for thermo-poro-mechanics (e.g., the Terzaghi, Mandel, and McNamee-Gibson problems), the opening of static fractures, and fracture propagations for viscosity and toughness dominated systems, matching the analytical solutions (e.g., Kim and Moridis (2013)).

Simulation domain

We take a 3D simulation domain, as shown in Fig.4. We discretize the domain for geomechanics with 50, 5, 50 gridblocks in x, y and z directions, respectively. The direction of the minimum compressive principal total stress is perpendicular to the x-z plane. The sizes of the gridblocks in the x and z direction are uniform, i.e., Δx = Δz = 3m. The sizes of the gridblocks in the y direction are non-uniform, i.e., 0.1m, 0.5m, 3.0m, 10.0m, and 20.0m. We take 12GPa of Young's modulus, E, and 0.3 of Poisson's ratio, v, respectively. The tensile strength of shale for the reference case is Tc = 10MPa. We take β = 1.0for contribution of effective shear stress to tensile failure.

Fig. 4

A domain of hydraulic fracturing simulation in 3D. We assume a simultaneous fracturing operation, where the horizontal wells are aligned with the direction of the minimum compressive principal total stress, Sh.

Fig. 4

A domain of hydraulic fracturing simulation in 3D. We assume a simultaneous fracturing operation, where the horizontal wells are aligned with the direction of the minimum compressive principal total stress, Sh.

Close modal

For flow, we share the same discretized domain of geomechanics except an additional layer in the y direction, having 50, 6, 50 gridblocks in x, y and z directions. The additional layer, the thickness of which is 1.0m, represents half of the width of the stimulated reservoir zone shown in Fig. 1. We assume water to be injected at (x=75m, z=-1440m), which is an initial fractured node. We take 8.645 ×10−18m2 of intrinsic permeability for the rock matrix, where 1 Darcy is 9.87 ×10−13m2. For relative permeability, we use Sr,w = 0 08 and Sr,g = 0.01. np = 3 0 and ac = 0.14 are taken for the modified cubic law. Biot's coefficient, α, is 1.0, and the thermal dilation coefficient is αT= 4.5×10−5°C−1. We apply the dual continuum approach to the simulated zone, which consists of the main fracture and the damaged zone that includes small local fractures. The permeability of the main fracture, kfp, is calculated from Eq. 3, and that of the damaged zone, kmp, is calculated from kmp=χmkfp, where χm = 10−3. The lower limit of kpf is 60mD. For Eq. 5, we assume Kl = 500MPa for the both continua. We take ηf = 0.01and ηm = 0.99

For the initial conditions, the reservoir pressure is initially 17.10MPa at 1350m in depth with the 12.44kPa/m gradient. Initial temperature is 58.75°C at 1350 m in depth with the 0.025°C /m gradient. Based on estimation of the generic total stress distribution (Zhang and Stephansson, 2010), the initial total principal stresses are −36.40MPa, and −23.30MPa, and −29.12MPa at 1350m in depth in x, y, and z directions, respectively, where tensile stress is positive. The corresponding stress gradients are −27.0kPa/m, −17.59kPa/m, and −21.57kPa/m, respectively. The injection rate is 40kg/s at the injection point. We have no-flow boundary condition for flow. For geomechanics, there are no-horizontal displacement boundary conditions at sides, except the fractured nodes, and have no displacement boundary at the bottom. We also consider gravity and the bulk density is 2200kg/m3.

For capillarity, we have the van Genutchen capillary pressure model (Genuchten, 1980), written as

Pc=Πc((SwSr,w1Sr,gSr,w)1λp1)1λp,
(6)

where pc, λp, and Πc are capillary pressure, the exponent that characterizes the capillary pressure curve, and the capillary modulus, respectively. We take λp = 0.45 and Πc = 20kPa in Eq. 6. Sr,w = 0.05 and Sr,g = 0.0 are taken for capillarity, which are slightly smaller than those used in the relative permeability model, in order to prevent unphysical behavior (Moridis et al., 2008).

Numerical results

Fig. 5 shows fracture propagation and fracture opening at the main fracture at different times. Due to the point source of water injection, the fracture propagates in an elliptic shape, as shown in the left figure. From the righ figure, the aperture of the fracture is continuous when the fracture is open. The fracture propagation is stable, which means that it depends on the injection time, and that it can be controlled by the injection rate and time.

Fig. 5

Left: growth of the fracture induced by water injection. Right: the corresponding fracture opening. The fracture propagation depends on the injection time.

Fig. 5

Left: growth of the fracture induced by water injection. Right: the corresponding fracture opening. The fracture propagation depends on the injection time.

Close modal

Fig. 6 shows flow and geomechanical responses at different monitoring points. Fracturing occurs more at early times than at late times. Due to discontinuity of fracturing events in time, we identify pressure fluctuation under the constant rate of water injection, observing that pressure drops every time when failure occurs. Fluctuation of pressure is much severer at early times than at late times, because the fracture volume is small at early times. Pressure responses of incompressible fluid are more sensitive to smaller fractured areas. Then the fluctuation decreases as the fracture length grows, because the fractured area becomes large. This pressure oscillation can also be found in Tzchichholz and Herrmann (1995). Along with the pressure fluctuation, the aperture of the fracture has the same behavior as pressure does. For displacement, we observe the uplift at the top of the simulation domain, because of expansion of the domain induced by fluid injection.

Fig. 6

Evolution of geomechanics and flow variables: (a) number of fractured nodes, (b) pressure at the injection point (b), (c) aperture of the fracture at the injection point, (d) displacement at (x=75m, z=-1350m). Fracturing events are discontinuous, which results in fluctuation in pressure and the aperture. The oscillation is more dominant at early times due to small volume of the stimulated zone.

Fig. 6

Evolution of geomechanics and flow variables: (a) number of fractured nodes, (b) pressure at the injection point (b), (c) aperture of the fracture at the injection point, (d) displacement at (x=75m, z=-1350m). Fracturing events are discontinuous, which results in fluctuation in pressure and the aperture. The oscillation is more dominant at early times due to small volume of the stimulated zone.

Close modal

In Fig. 7, we find existence of reservoir gas within the fracture. In particular, we identify a large amount of gas within the fracture at late times. As the fracture grows, reservoir gas initially captured in the stimulated zone moves upward due to buoyancy, taking the upper part of the fracture, while the injected water fills the lower part of the fracture because of higher density. In Fig. 8, we also find that, at early times, the water saturation front is almost identical to the fracture tip, showing that the fracture is mostly filled with injected water. However, at late times, advance of the water front is retarded, compared to the fracture propagation, yielding a significant gap between the water front and the fracture top, which is filled with reservoir gas. The fracture propagation is not the same as propagation of the water front, because fracturing is governed by geomechanics, whereas water saturation is determined by fluid flow, just as anticipated in the previous section. We thus cannot assume the fracture to be fully saturated with injected water when estimating dimension of the fracture.

Fig. 7

Gas saturation at different times. The color-bar indicates gas saturation. Gas is accumulated at the upper part of the fracture, while water fills the lower part.

Fig. 7

Gas saturation at different times. The color-bar indicates gas saturation. Gas is accumulated at the upper part of the fracture, while water fills the lower part.

Close modal
Fig. 8

Distributions of water saturation, aperture, dimensionless pressure (Pd) at x=75m.pd=ppLpUpL, where pu=40MPa and PL=17.1MPa. The aperture is continuous along the fracture, while gas saturation is somewhat discontinuous. At early times, water fills almost all part of the fracture, whereas gas takes considerable volume of the fracture at late times.

Fig. 8

Distributions of water saturation, aperture, dimensionless pressure (Pd) at x=75m.pd=ppLpUpL, where pu=40MPa and PL=17.1MPa. The aperture is continuous along the fracture, while gas saturation is somewhat discontinuous. At early times, water fills almost all part of the fracture, whereas gas takes considerable volume of the fracture at late times.

Close modal

In Fig. 9, we find considerable leak-off of water to the reservoir. Gas saturations at the damaged zone and intact reservoir nearest to the stimulated zone are significantly low, indicating that considerable amount of water infiltrates into the reservoir. When the leak-off water cannot be removed, it might potentially reduce productivity of gas production.

Fig. 9

Distribution of gas saturations at the damaged zone (left figure) and the rock matrix (right figure). The large amount of water infiltrates from the fracture into the reservoir.

Fig. 9

Distribution of gas saturations at the damaged zone (left figure) and the rock matrix (right figure). The large amount of water infiltrates from the fracture into the reservoir.

Close modal

Finite element solution to Maxwell's equations for electromagnetic modeling

EM geophysical responses to an Earth model can be simulated by solving Maxwell's equations. Maxwell's equations in the frequency domain are given by

×e(r)=i^ωb(r),
(7)
μo1×b(r)=σ(r)e(r)+i^εoωe(r)+Js(r),
(8)
e(r)=ρvεo1,
(9)
b(r)=0,
(10)

where e(r) is the electric field at position r, Js(r)is an electric current source at angular frequency ω, μo is the magnetic permeability, σ(r) is the electrical conductivity at position r, εo is the dielectric permittivity, and ρv is the electric charge density. In the paper, μo within the earth and air is assumed to be constant and is set to that of free space (4π × 10−7 H/m). εo is generally a tensor quantity but assumed to be a constant scalar quantity (8.85×10−12 F/m). σ(r) is also assumed to be a scalar quantity. i^ indicates the imaginary unit.

By combining Eq. 7 with Eq. 8, we have the electric-field full-wave equation,

1μo××e(r)εoω2e(r)+i^ωσe(r)+i^ωJs(r)=0.
(11)

Eq. 11 is modified to the electric-field diffusion equation by neglecting its second term. By doing this, we have:

1μo××e(r)+i^ωσe(r)+i^ωJs(r)=0.
(12)

We consider the finite-element (FE) formulation for solving the electric-field diffusion equation. The FE formulation of Eq. 12 starts with its equivalent weak statement (Um et al., 2013). The development of the weak statement requires the multiplication of Eq. 6 by the edge basis function (Nédélec, 1980) and the integration over the model domain of V, resulting in

Venie(r)(1μo××ee(r)+i^ωσeee(r)+i^ωJs(r)dV=0,
(13)

where we use the first-order edge basis functions and tetrahedral elements. The superscript e denotes the eth tetrahedral element, nie(r) with i varying from 1 to 6 is a set of edge basis functions and Ve is the volume of the eth tetrahedral element. The set of nie(r) used in Eq. 7 is also chosen as the basis set. Thus, the electric field at a point inside or on a given element is expanded as

ek(r)=j=16ejk(r)=j=16ujknjk(r),
(14)

where ujk is the unknown amplitude of the electric field on edge j of the kth element.

Substituting Eq. 14 into Eq. 13 and using vector calculus identities with the homogeneous Dirichlet boundary conditions yield

(Ae+μoωi^Be)ue=μωi^se,
(15)

Where

the (i,j)th  element of Ae=Ve×nie(r)×nje(r)dV,
(16)
the (i,j)th  element of Be=Venie(r)σenje(r)dV
(17)
the ith  element of se=Venie(r)Js(r)dV,
(18)
ue=[u1eu2eu6e]
(19)

Eq. 15 is considered local, because it results from integration over each individual element. Based on connectivity information about elements in V, the local systems of FE equations derived from the individual elements are assembled into a single global system of equations. Then, the superscript e is dropped. The solution of the global version of Eq. 15 determines the electric fields in V. The magnetic fields are interpolated from the electric fields using the Faraday's law. The system matrix of Eq. 15 is unstructured, sparse, complex and symmetric. Eq. 15 can be solved using an LU decomposition (Golub and Van Loan, 1996) or quasi-minimal residual method (Freund, 1992).

Electrical conductivity models

Once a set of reservoir parameters (e.g., porosity ϕ, pore fluid saturation Sw, and electrical conductivity of pore fluid σw) are determined from coupled flow-geomechanical simulation, the buck electrical conductivity (σ) of a shaly reservoir can be estimated by Waxman-Smits equation (Waxman and Smits, 1968):

σ=ϕmSwn[σw+BQvSw1],
(20)
B=4.78×108[10.6eσw0.0013],
(21)
Qv=ρg(1ϕ)ϕ1CEC.
(22)

In this study, the rock physics parameters of Eqs. 20, 21 and 22 are empirically chosen; both the cementation (m) factor and the saturation (n) factor are set to 2; the grain density (ρg) is set to 2,370 (kg/m3); the cation exchange capacity (CEC) is set to 38,000 (Coulomb/kg); the electrical conductivity of the pore fluid is set to 3.33 (S/m). Thus, Eq. 20 transforms the reservoir parameters to an electrical conductivity value.

Hydraulic fracturing operations force fluids into the formation under high pressure, producing highly-localized changes in porosity and pore fluids. Because EM geophysical methods are sensitive to porosity, fluid saturation and chemistry of the fluid in the pore space, they are a promising tool for illuminating migration pathways of the injected fluids and can complement microseismic methods.

A primary goal of our EM modeling analysis is to examine capabilities of an EM geophysical method for sensing bulk conductivity changes that are directly related to migration pathways of injected fluids. For successful EM detection for the fluid movement, it is important to ensure a sufficient contrast in electrical conductivity between intact and the stimulated zones. The sufficient contrast can be realized by injecting the brine (3.3 S/m). One can also consider electrically-engineered injection fluids that have high conductivity values (e.g. 10,000 S/m). The fabrication of a mixture of such fluids and proppant is currently feasible. For detailed information about properties of electrically-engineered fluids and proppant, a reader is referred to Natália Saliés (2012) and references therein. Magnetically-engineered fluids are also feasible. One can find more information about them from Borglin et al. (2000) and Oldenburg et al. (2000).

Based on Waxman-Smits equation, we transform reservoir parameters to realistic electrical conductivity models. Fig. 10 shows cross-sectional views of reservoir parameters (the first column for reservoir porosity and the second column for water saturation) and electrical conductivity models (the third column) after hydraulic fracturing starts along with a high-conductivity fluid (10,000 S/m). In Fig. 10, the porosity and the saturation are volume-averaged on the stimulated zone. For a comparison purpose, we also generate the electrical conductivity models (the fourth column) when the brine (3.33S/m) is injected. As the hydraulic fracturing operation continues, perturbations in porosity and saturation gradually increase in all the directions (e.g., leak-off of water).

Fig. 10

Cross-sectional view (x-z plane at the first y layer, stimulated zone) of reservoir models and resultant electrical conductivity models. The 1st column and the 2nd column show changes in porosity and fluid saturation values over time, respectively. The 3rd column shows changes in conductivity values when the electrically-engineered high-conductivity fluid (10,000 S/m) is injected. The 4th column shows changes in conductivity values when the brine (3.33 S/m) is injected. At the top of the third column, the yellow vertical line segments indicate observation wells for the vertical crosswell EM configuration. Two yellow ⊗ indicates the positions of horizontal wells for the horizontal crosswell EM configuration. The color charts for the 1st and 2nd column are linear and those for the 3rd and 4th column are based on the common log.

Fig. 10

Cross-sectional view (x-z plane at the first y layer, stimulated zone) of reservoir models and resultant electrical conductivity models. The 1st column and the 2nd column show changes in porosity and fluid saturation values over time, respectively. The 3rd column shows changes in conductivity values when the electrically-engineered high-conductivity fluid (10,000 S/m) is injected. The 4th column shows changes in conductivity values when the brine (3.33 S/m) is injected. At the top of the third column, the yellow vertical line segments indicate observation wells for the vertical crosswell EM configuration. Two yellow ⊗ indicates the positions of horizontal wells for the horizontal crosswell EM configuration. The color charts for the 1st and 2nd column are linear and those for the 3rd and 4th column are based on the common log.

Close modal

It is noteworthy that because the fracture propagation is not the same as the injected fluid movement, the anomalous conductivity distribution gradually deviates over time from the anomalous porosity distribution and increasingly resembles the anomalous saturation distribution. Accordingly, EM responses to the given conductivity models will be more directly correlated with the migration pathways of the injected fluids but might be marginally sensitive to unsaturated fractures. Therefore, when EM geophysical methods are used to characterize hydraulically-induced fractures, one should be aware of the possibility that the unsaturated fractures may not be imaged. In contrast, as mentioned earlier, the microseismic methods are based on seismic events associated with fracturing but are insensitive to saturation changes in late times. Therefore, we can expect that the joint analysis of microseismic and EM data would significantly reduce the ambiguity in characterizing hydraulically-induced fractures. Also note that the use of the high-conductivity fluid does not change the overall geometry of anomalous conductivity distribution (compare the third column of Fig. 10 with the fourth column) but significantly increase the conductivity contrast over time. Thus, as will be demonstrated in the next section, its use makes EM geophysical methods more suitable for monitoring migration pathways of the injected fluids.

Crosswell EM methods

To monitor migration pathways of the injected fluids, we consider crosswell EM methods. The crosswell EM methods (Fig. 11) interrogate electrical conductivity structures between wells and can yield detailed 2D/3D images. The EM methods typically employ solenoids as a transmitter (a wire coil carrying alternating currents). The emitted magnetic source fields are called primary magnetic fields (Bo) and induce currents in nearby formations. In turn, the induced currents generate secondary magnetic fields. The primary and secondary magnetic fields are measured in the other well with magnetic receivers. After a number of sources are excited at different positions in one well and the magnetic fields are measured at a number of receiver positions in the other well, the conductivity structures between the wells are imaged through inverse modeling. For details of the crosswell EM methods and their applications, a reader is referred to Alumbaugh and Morrison (1995), Wilt et al. (1995), Zeng et al. (2000), and Gao et al. (2008). To successfully sense and monitor migration pathways of the injected fluids, conductivity changes due to their migration should produce measurable perturbation in the magnetic fields. In this modeling study, we focus on the sensitivity analysis of the crosswell EM methods for the conductivity changes.

Fig. 11

A typical crosswell EM survey configuration (Schlumberger, 2006). Sources and receivers are lowered into the wells. The primary magnetic fields of the transmitter coil induce the electrical current (yellow). The current produces the secondary magnetic fields (pink). The receivers record the both primary and magnetic fields.

Fig. 11

A typical crosswell EM survey configuration (Schlumberger, 2006). Sources and receivers are lowered into the wells. The primary magnetic fields of the transmitter coil induce the electrical current (yellow). The current produces the secondary magnetic fields (pink). The receivers record the both primary and magnetic fields.

Close modal

We consider vertical and horizontal crosswell EM configurations. The proposed configurations are shown in Fig. 10a. The vertical crosswell EM configuration intends to better sense the height and the width of fractures between wells. The horizontal crosswell EM configuration would be better sensitive to the thickness of the fractures. For both crosswell EM configurations, well spacing is set to 100 m. A source frequency is set to 3,000 Hz. To allow the high frequency source in a well environment, it is assumed that a part of casing where the sources and the receivers are installed is cased with non-metallic materials (e.g. fiberglass). It is also assumed that magnetic receivers can detect 1% difference in the magnetic field amplitude and the receiver noise level is 1.0−10 (Amp/m).

The proposed crosswell EM configurations have been simulated using a 3D FE electromagnetic simulator (Um et al., 2013) over the electrical conductivity models shown in Fig. 10. Vertical crosswell EM responses to the conductivity models with the high-conductivity fluid and the brine are shown at selected source positions in Figs. 12 and 13, respectively. Horizontal crosswell EM responses to the conductivity models with the high-conductivity fluid and the brine are shown in Figs. 14 and 15, respectively. For simplicity, the total magnetic fields and their relative difference with respect to the magnetic fields of 0.02S/m half-space (i.e. the background conductivity before the injection) are plotted. First, note that all magnetic fields (the left column of Figs. 12 through 15) are above the noise level prescribed above. As the fracture grows larger, both horizontal and vertical crosswell measurements start to sense the fracture. When the brine is injected, the perturbation in the magnetic field is marginal (Figs. 13 and 15). However, this does not imply that the brine is an ineffective tracer when the crosswell EM method is employed to map saturated fractures. In fact, the brine movements have been successfully detected in many enhanced oil recovery operations. However, in this particular case, the brine injection produces small changes in the conductivity (the fourth column of Fig. 10) due to limited increase of porosity from 0.049 to 0.059 during hydraulic fracturing, resulting in the marginal perturbation in the magnetic field. In contrast, Figs. 12 and 14 show that the use of the high-conductivity fluid significantly (about an order of magnitude) improves the sensitivity of the crosswell EM measurements to conductivity changes due to hydraulic fracturing operations. We conclude that along with an electrically-engineered fluid, the crosswell EM method can serve as an effective mapping tool for conductivity changes that are directly related with the migration pathways of the injected fluid.

Fig. 12

Magnetic field amplitudes (the left column) and their relative difference (the right column) with respect to the background magnetic field when the vertical crosswell configuration is applied to the conductivity model shown in the third column of Fig. 10. Receivers are placed along the axis that is parallel to the z-axis and passes through the point (x=130 m and y=0 m).

Fig. 12

Magnetic field amplitudes (the left column) and their relative difference (the right column) with respect to the background magnetic field when the vertical crosswell configuration is applied to the conductivity model shown in the third column of Fig. 10. Receivers are placed along the axis that is parallel to the z-axis and passes through the point (x=130 m and y=0 m).

Close modal
Fig. 13

Magnetic field amplitudes (the left column) and their relative difference (the right column) with respect to the background magnetic field when the vertical crosswell configuration is applied to the conductivity model shown in the fourth column of Fig. 10. Receivers are placed along the axis that is parallel to the z-axis and passes through the point (x=130 m and y=0 m).

Fig. 13

Magnetic field amplitudes (the left column) and their relative difference (the right column) with respect to the background magnetic field when the vertical crosswell configuration is applied to the conductivity model shown in the fourth column of Fig. 10. Receivers are placed along the axis that is parallel to the z-axis and passes through the point (x=130 m and y=0 m).

Close modal
Fig. 14

Magnetic field amplitudes (the left column) and their relative difference (the right column) with respect to the background magnetic field when the horizontal crosswell configuration is applied to the conductivity model shown in the third column of Fig. 10. Receivers are placed along the axis that is parallel to the y-axis and passes through the point (x=130 m and z=−1440 m).

Fig. 14

Magnetic field amplitudes (the left column) and their relative difference (the right column) with respect to the background magnetic field when the horizontal crosswell configuration is applied to the conductivity model shown in the third column of Fig. 10. Receivers are placed along the axis that is parallel to the y-axis and passes through the point (x=130 m and z=−1440 m).

Close modal
Fig. 15

Magnetic field amplitudes (the left column) and their relative difference (the right column) with respect to the background magnetic field when the horizontal crosswell configuration is applied to the conductivity model shown in the fourth column of Fig. 10. Receivers are placed along the axis that is parallel to the y-axis and passes through the point (x=130 m and z=−1440 m).

Fig. 15

Magnetic field amplitudes (the left column) and their relative difference (the right column) with respect to the background magnetic field when the horizontal crosswell configuration is applied to the conductivity model shown in the fourth column of Fig. 10. Receivers are placed along the axis that is parallel to the y-axis and passes through the point (x=130 m and z=−1440 m).

Close modal

We investigated geomechanical and flow responses during hydraulic fracturing operations induced by water-injection. From the numerical results, we identified that fracture propagation is not the same as propagation of the water front. The gap between the water front and the fracture top is more dominant at later times when a large amount of water is injected, and the gap is filled with reservoir gas. We also found considerable leak-off of water to the reservoir. The inconsistency between the fracture volume and the volume of injected water might underestimate the fracture propagation, when we take the simple assumption that the fracture is fully saturated with injected water. We identified fluctuation in pressure and displacement under constant water injection. As the fracture dimension becomes larger, the fluctuation decreases. All the complex physical processes above are fundamentally based on tight coupling between flow and geomechanics, which have different time scales between them. Thus, rigorous modeling of coupled flow and geomechanics is strongly recommended to capture accurate reservoir behavior.

Owing to rigorous 3D modeling capabilities of the recently-developed coupled flow-geomechanical simulator, T+M, we have been able to generate a series of realistic electrical conductivity models that reveals transient changes of conductivity in shaly reservoirs during hydraulic fracturing operations. Because the injected fluid movement is separated from the fracture propagation in late times, the conductivity models are more closely related with the migration pathways of the injected fluids rather than unsaturated fractures. The 3D FE Maxwell's equations solver was used in this study to simulate crosswell EM methods over a series of the conductivity models. The modeling analysis suggests the crosswell EM method can serve as an effective mapping tool for conductivity changes along the migration pathways of the injected fluids. Since our sensitivity analysis demonstrates sufficient perturbation in the magnetic fields during hydraulic fracturing, we expect that 3D inverse modeling of the crosswell EM measurements will determine spatial distribution of electrical conductivity around the wells and make it possible to delineate fluid migration pathways and fracture geometry. The resulting conductivity images will also provide a clearer understanding of the fractured reservoir. Further sensitivity analysis with different types of electromagnetically-engineered fluids and proppants, optimal geophysical survey designs for typical hydraulic fracturing scenarios and subsequent geophysical imaging experiments will be performed in our future research.

This paper was selected for presentation by an SPE program committee following review of information contained in an abstract submitted by the author(s). Contents of the paper have not been reviewed by the Society of Petroleum Engineers and are subject to correction by the author(s). The material does not necessarily reflect any position of the Society of Petroleum Engineers, its officers, or members.

This study was supported by the US Environmental Protection Agency, Office of Water, under an Interagency Agreement with the U.S. Department of Energy at the Lawrence Berkeley National Laboratory through Contract No. DE-AC02-05CH11231, and by RPSEA (Contract No. 08122-45) through the Ultra-Deepwater and Unconventional Natural Gas and Other Petroleum Resources Research and Development Program as authorized by the US Energy Policy Act (EPAct) of 2005. The research described in this article has been funded wholly (or in part) by the U.S. Environmental Protection Agency through Interagency Agreement (DW-89-92235901-C) to the Lawrence Berkeley National Laboratory. The views expressed in this article are those of the author(s) and do not necessarily reflect the views or policies of the EPA.

Adachi
,
J.
,
Siebrits
,
E.
,
Peirce
,
A.
,
Desroches
,
J.
,
2007
.
Computer simulation of hydraulic fractures.
Int. J. Rock. Mech. Min.
,
44
,
739
757
.
Alumbaugh
,
D.
,
Morrison
,
H.F.
1995
.
Monitoring subsurface changes over time with cross well electromagnetic tomography.
Geophysical Prospecting
43
:
873
902
.
Berryman
JG.
2002
Extension of poroelastic analysis to double-porosity materials: New technique in microgeomechanics.
J. Eng. Mech. ASCE
;
128
(
8
):
840
847
.
Borglin
,
S.E.
,
G.J.
Moridis
,
C.M.
Oldenburg
,
2000
,
Experimental Studies of the Flow of Ferrofluids in Porous Media
,
Transport in Porous Media,
41
:
61
80
.
Cipolla
,
C.L.
,
Lolon
,
E.P.
,
Erdle
,
J.C.
,
Rubin
.
B.
,
2010
Reservoir modeling in shale-gas reservoirs.
SPE Reserv. Eval. Eng.
638
653
.
Coussy
,
O.
,
2004
Poromechanics.
Chichester, England
:
John Wiley and Sons.
Freund
,
R.
,
1992
,
Conjugate gradient type methods for linear systems with complex symmetric coefficient matrices
,
SIAM Journal Scientific Statistical Computing,
13
:
425
448
.
Gao
,
G.
,
D.
Alumbaugh
,
P.
Zhang
,
H.
Zhang
,
C.
Levesque
,
R.
Rosthal
,
J.
Liu
,
A.
Abubakar
, and
T.
Habashy
,
2008
,
Practical Implications of Nonlinear Inversion For Cross-well Electromagnetic Data Collected In Cased-wells
,
SEG Annual Meeting Abstract
.
Golub
,
G. H.
and
C.F.
Van Loan
,
1996
,
Matrix computation
, 3rd edition:
Johns Hopkins University Press
.
Ji
,
L.
, and
Settari
,
A.
, and
Sullivan
,
R.B.
2009
A novel hydraulic fracturing model fully coupled with geomechanics and reservoir simulation.
SPEJ..
423
430
.
Johnson
,
R.
and
J.
Shrallow
,
2011
,
Ambiguity in microseismic monitoring
,
SEG Annual Meeting Abstract
.
Dean
,
R.H.
and
Schmidt
,
J.H.
,
2008
Hydraulic fracture predictions with a fully coupled geomechanical reservoir simulation.
SPEJ,
14
(
4
):
707
714
Fisher
,
K.
and
Warpinski
N
2012
Hydraulic fracture-height growth: real data.
SPE Prod. Oper.
27
(
1
):
8
19
Fjaer
,
E.
,
Holt
,
R.M.
,
Horsrud
,
P.
,
Raaen
,
A.M.
,
Risnes
,
R.
,
2008
.
Petroleum related rock mechanics,
2nd ed.
Elsevier
,
B.V., Amsterdam, The Netherlands
.
Kargbo
,
D.M.
,
Wilhelm
,
R.G.
,
Campbell
,
D.J.
,
2010
Natural gas plays in the Marcellus shale: Challenges and potential opportunities.
Environ. Sci. Technol
.,
44
(
5
):
5679
5684
.
Kim
,
J.
,
Sonnenthal
,
E.
, and
Rutqvist
,
J.
2012
Formulation and sequential numerical algorithms of coupled fluid/heat flow and geomechanics for multiple porosity materials.
Int. J. Numer. Meth. Engrg.
92
:
425
456
, doi: .
Kim
J
, and
Moridis
,
G.J.
2013
Development of the T+M coupled flow-geomechanical simulator to describe fracture propagation and coupled flow-thermal-geomechanical processes in tight/shale gas systems.
Comput. Geosci.
60
:
184
198
.
King
,
G. E.
,
2012
Hydraulic fracturing 101: What every representative, environmentalist, regulator, reporter, investor, university researcher, neighbor and engineer should know about estimating frac risk and improving frac performance in unconventional gas and oil wells.
SPE Hydraulic Fracturing Technology Conference
,
The woodland, TX
,
Feb. 6-8. 2012
.
Moridis
,
G.J.
,
Kowalsky
,
M.B.
,
Pruess
,
K.
,
2008
Tough+hydrate v1.0 user's manual: A code for the simulation of system behavior in hydrate-bearing geologic media.
In: Report LBNL-00149E.
Lawrence Berkeley National Laboratory
,
Berkeley, CA
.
Nédélec
,
J.-C.
,
1980
,
Mixed finite elements in R3
:
Numerische Mathematik,
35
:
315
341
.
Nordren
,
R.P.
,
1972
.
Propagation of a vertical hydraulic fracture.
SPEJ.
12
(
8
),
306
314
,
SPE7834
.
Oldenburg
,
C.
,
S.E.
Borglin
, and
G. J.
Moridis
,
2000
,
Numerical Simulation of Ferrofluid Flow for Subsurface Environmental Applications
,
Transport in Porous Media,
38
:
319
344
.
Osborn
,
S.G.
,
Vengosh
,
A.
,
Warner
,
N.R.
,
Jackson
,
R.B.
,
2011
,
Methane contamination of drinking water accompanying gas-well drilling and hydraulic fracturing.
PNAS
,
108
:
8172
8176
.
Page
,
J.C.
,
Miskimins
,
J.L.
,
2009
,
A comparison of hydraulic and propellant fracture propagation in a shale gas reservoir.
J. Can. Petrol. Technol.
,
48
(
5
):
26
30
.
Perkins
,
T.K.
,
Kern
,
L.R.
,
1961
,
Widths of hydraulic fractures.
JPT
13
(
9
),
937
949
,
SPE89
.
Ribeiro
,
L. H.
,
Sarma
,
M. M.
,
2013
Fluid selection for energized fracture treatments.
SPE Hydraulic Fracturing Technology Conference
,
The woodland, TX
,
Feb. 4-6. 2013
.
Ruiz
,
G.
,
Ortiz
,
M.
,
Pandolfi
,
A.
,
2000
,
Three-dimensional finite-element simulation of the dynamic Brazilian tests on concrete cylinders.
Int. J. Numer. Meth. Engrg.
48
,
963
994
.
Rutqvist
,
J.
,
Rinaldi
,
A.P.
,
Cappa
,
F.
,
Moridis
,
G.J.
,
2013
,
Modeling of fault reactivation and induced seismicity during hydraulic fracturing of shale-gas reservoirs.
,
J. Pet. Sci. Eng.
107
:
31
44
Rutqvist
,
J.
,
Stephansson
,
O.
.
2003
,
The role of hydromechanical coupling in fractured rock engineering.
Hydrogeol. J,
.
11
:
7
40
.
Saliés
,
Natália
Gastão
,
2012
,
Study on the feasibility of using electromagnetic methods for fracture diagnostics
,
Thesis of Master of Science in Engineering
,
University of Texas at Austin
.
Schlumberger
,
2006
,
Middle East and Asia Reservoir Review
,
24
32
Snow
,
D.T.
,
1965
,
A parallel plate model of fractured permeable media.
,
Ph. D thesis
,
University of California
,
Berkeley
.
Tzchichholz
,
F.
,
Herrmann
,
H. J.
,
1995
,
Simulations of pressure fluctuations and acoustic emission in hydraulic fracturing.
Phys. Rev. E,
51
(
3
):
1961
1970
.
Um
,
E.
,
M.
Commer
,
G.
Newman
,
2013
,
Efficient pre-conditioned iterative solution strategies for the electromagnetic diffusion in the Earth: finite-element frequency-domain approach
,
Geophysical Journal International,
193
:
1460
1473
.
van
Genuchten
,
1980
,
A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.
Soil Sci. Soc. Am, J.
, (
44
):
892
898
.
Vermylen
,
J.P
and
Zoback
,
M.D.
2011
Hydraulic fracturing, microseismic magnitudes, and stress evolution in the Barnett Shale.
USA. Hydr. Frac. Tech. Conf.
The woodland, TX
,
24–26 Jan
.
Warpinski
,
N. R.
,
Kramm
,
R. C.
,
Heinze
,
J. R.
,
Waltman
,
C. K.
,
2005
,
Comparison of single- and dual-array microseismic mapping techniques in the Barnett shale.
SPE ATCE
,
Dallas
,
Oct
.
9
12
.
2005. SPE95568
Waxman
,
M. H.
and
Smits
,
L. J. M.
,
1968
,
Electrical conductivities in oil-bearing shaly sands
,
SPEJ
8
:
107
122
.
Wilt
,
M.
,
Alumbaugh
,
D.
,
Morrison
,
H.F.
,
Becker
,
A.
,
Lee
,
K.H.
,
Deszcz-Pan
,
M.
1995
,
Crosshole electromagnetic tomography: System design considerations and field results.
Geophysics
60
:
871
885
.
Witherspoon
,
P.A.
,
Wang
,
J.S.Y.
,
Iwai
,
K.
,
Gale
,
J.E.
,
1980
,
Validity of cubic law for fluid flow in a deformable rock fracture.
Water. Resour. Res.
16
(
6
):
1016
1024
.
Zang
,
A.
,
Stephansson
,
O.
,
2010
,
Stress field of the Earth's crust.
Springer Science+Business Media, B.V.
Zeng
,
W.
,
Zhao
,
W.
,
Wilt
,
M.
2000
.
Three Successful Large-scale Field Tests of Crosshole Electromagnetic Tomography Technology.
Paper SPE 63143
.
Zoback
,
M.
,
Kitasei
,
S.
,
Copithorne
,
B.
.
2010
Addressing the environmental risks from shale gas development.
Worldwatch Institute Briefing Paper 1
(
Worldwatch Inst
,
Washington DC
).
Zoback
,
M.D.
2007
Reservoir geomechanics.
Cambridge
,
Cambridge university press
.

Supplementary data