## Abstract

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.

## Introduction

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.

## Hydraulic Fracturing Simulation

### 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.

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^{−4}Pa·s and 1.78× 10^{−5}Pa·s, respectively. Density of gas at the reservoir condition, 1.15 ×10^{2}kg·m^{−3}, cannot also be neglected, compared with density of water, 9.89 × 10^{2} 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, *f _{w}*, after ignoring the gradient of capillary pressure within a fracture.

where *N _{g}, k_{r,J}*, and

*μ*are the gravity number, relative permeability and viscosity of phase

_{J}*J*. The subscript

*g*and

*w*indicate gas and water phases, respectively.

*S*and

_{J}*S*are saturation and residual (irreducible) saturation of phase

_{r,J}*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.

### 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

where *σ′ _{c}* is the effective stress for tensile failure, which consists of shear and normal effective stresses.

*t*′

*′*

_{n}, t*and*

_{t}*t*′

*are normal and shear effective stresses acting on a fracture plane (Fig. 3), and*

_{s}*T*is tensile strength. Using

_{c}*β*of Eq. 2, we account for contribution of shear effective stress to tensile failure.

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

where *ω _{f}* is the aperture, also called fracture opening.

*Q*is flow rate of water.

_{w}*H*is the fracture plate width, and

**g**is the gravity vector. For

*n*= 3.0, Eq. 3 reduces to the original cubic law.

_{p}*a*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.

_{c}### 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

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. *K _{l}, K*

_{s},

*η*, are the drained bulk modulus, intrinsic solid grain modulus, the volume fraction, the thermal dilation coefficient, and Biot's coefficient for subelement

_{l}, α_{T,l}α_{l}*l*within a gridblock, respectively.

*σ*is total mean (volumetric) stress at the gridblock.

_{v}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

where the superscript *n* is the time level in time discretization. *c* and *c _{T}* correspond to the pore-compressibility and thermal expansivity used in reservoir simulation, respectively. $\Delta \Phi 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* = 3*m*. 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 *T _{c}* = 10

*MPa*. We take

*β*= 1.0for contribution of effective shear stress to tensile failure.

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^{−18}*m*^{2} of intrinsic permeability for the rock matrix, where 1 Darcy is 9.87 ×10^{−13}*m*^{2}. For relative permeability, we use *S _{r,w}* = 0 08 and

*S*= 0.01.

_{r,g}*n*= 3 0 and

_{p}*a*= 0.14 are taken for the modified cubic law. Biot's coefficient,

_{c}*α*, is 1.0, and the thermal dilation coefficient is

*α*= 4.5×10

_{T}^{−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=\chi mkfp$, where

*χ*= 10

^{m}^{−3}. The lower limit of $kpf$ is 60mD. For Eq. 5, we assume

*K*= 500

_{l}*MPa*for the both continua. We take

*η*= 0.01and

_{f}*η*= 0.99

_{m}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/m^{3}.

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

where *p _{c}, λ_{p}*, and Π

*are capillary pressure, the exponent that characterizes the capillary pressure curve, and the capillary modulus, respectively. We take*

_{c}*λ*= 0.45 and Π

_{p}*= 20*

_{c}*kPa*in Eq. 6.

*S*= 0.05 and

_{r,w}*S*= 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).

_{r,g}### 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. 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.

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.

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.

## Electromagnetic Monitoring of Fluid Migration Pathways

### 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

where **e**(**r**) is the electric field at position **r, J**_{s}(**r**)is an electric current source at angular frequency *ω, μ _{o}* is the magnetic permeability, σ

**(r)**is the electrical conductivity at position

**r**,

*ε*is the dielectric permittivity, and

_{o}*ρ*is the electric charge density. In the paper,

_{v}*μ*within the earth and air is assumed to be constant and is set to that of free space (4π × 10

_{o}^{−7}H/m).

*ε*is generally a tensor quantity but assumed to be a constant scalar quantity (8.85×10

_{o}^{−12}F/m). σ(

**r**) is also assumed to be a scalar quantity. $i^$ indicates the imaginary unit.

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

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

where we use the first-order edge basis functions and tetrahedral elements. The superscript *e* denotes the e^{th} tetrahedral element, $nie(r)$ with *i* varying from 1 to 6 is a set of edge basis functions and *V ^{e}* is the volume of the e

^{th}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

where $ujk$ is the unknown amplitude of the electric field on edge *j* of the *k*^{th} element.

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

Where

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 S_{w}, 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):

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/m^{3}); 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).

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 (B_{o}) 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.

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.

## Conclusions

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.

## Acknowledgement

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.