Analytical solutions for flow to a well through a fault

Advances in Water Resources 29 (2006) 1790–1803 www.elsevier.com/locate/advwatres

Analytical solutions for flow to a well through a fault Erik I. Anderson

*

Department of Civil and Environmental Engineering, University of South Carolina, 300 Main Street, Columbia, SC 29208, United States Received 23 August 2005; received in revised form 20 December 2005; accepted 20 December 2005 Available online 14 February 2006

Abstract We present explicit analytical solutions to problems of steady groundwater flow to a pumping well in an aquifer divided by an infinite, linear fault. The transmissivity of the aquifer is allowed to jump from one side of the fault to the other to model the juxtaposition of host rocks with different hydrologic properties caused by faulting. The fault itself is represented as a thin anisotropic inhomogeneity; this allows the fault to act as a combined conduit–barrier to groundwater flow, as is commonly described in the literature. We show that the properties of the fault may be represented exactly by two lumped parameters—fault resistance and fault conductance—and that the effects of the fault on flow in the adjacent aquifer is independent of the fault width. We consider the limiting cases of a purely leaky and a purely conductive fault where the fault domain may be replaced exactly by internal boundary conditions, and we investigate the effects of fault properties on the flow behavior in the adjacent aquifers. We demonstrate that inferring fault properties based on field observations of head in the aquifer is inherently difficult, even when the fault may be described by one of the two limiting cases. In particular, the effects of a leaky fault and a conductive fault on heads and discharges in the aquifer opposite the fault from the well, are shown to be identical in some cases. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Groundwater; Fluid flow; Analytic; Leaky boundary; Faults

1. Introduction Faults and fault zones impact fluid flow by altering the local permeability of an aquifer and by placing rocks of different properties on either side of the fault. The study of fluid movement through fault zones has environmental and economic applications including the evaluation of water supply and quality [26,22,33,5], the evaluation of waste repositories [2,30,13], petroleum exploration [4,21, 32,24], and gold exploration [45]. Of particular environmental concern, faults can result in rapid transfer of contaminants between vertically adjacent aquifers or allow for rapid horizontal movement within the fault itself [12,10]. Numerous field studies indicate that faults may act as barriers to normal flow, conduits to flow within the fault plane, or a combination of both. However, evalu-

*

Fax: +1 803 777 0670. E-mail address: [email protected]

0309-1708/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2005.12.010

ating the bulk hydrologic properties of a fault zone through field observation and testing, and numerical simulation has proven difficult. Currently there is consensus on a simple qualitative model of the combined architecture and hydrologic properties of faults formed in the brittle, upper-crust; the model is applicable on a wide variety of scales and geologic environments. Caine et al. [14] described the architecture of a fault as consisting of three basic elements: first, the low permeability fault core is the location of most displacement and is characterized by the presence of breccia and/or cataclasite and gouge; second, the surrounding damage zone is characterized by an enhanced permeability due to dense fracturing; third, the protolith or host rock has permeabilities unaltered by the faulting, but may have different properties on either side of the damage zone due to the shifting of host rock when the faulting occurred. This layered feature of a low permeability fault zone layer lying between two high permeability layers of damage zone is a discrete model of an anisotropic medium.

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

Sigda et al. [41] identify two end-member fault types based on this architecture. In the first, the core zone is absent while the damage zone is well developed, and the fault acts solely as a conduit for fluid flow within the fault plane; we refer to this fault type here as a conductive fault. In the second, the damage zone is absent but the core is well developed and the fault acts solely as a barrier to fluid flow normal to the fault; we refer to this fault type as a leaky fault. Varying combinations of core and damage zones result in a fault that acts as a combined conduit–barrier system; we refer to this fault type as a general fault. The qualitative description of Caine et al. [14] is consistent with the results of quantitative studies performed on faults in granite [20,18], sandstone [38,4,21,43], and unconsolidated sediments [9,10]. Evans et al. [18] performed permeability tests on granitic rocks cored from a fault zone at different orientations, and report different properties in the protolith, the damaged zone, and the fault core. Their conceptual model of the fault system is of a narrow zone of anisotropic permeability with low permeability in the principal direction normal to the fault plane, high permeability in the principal direction parallel to the slip vector, and intermediate permeability in the fault plane normal to the slip vector. Antonellini and Aydin [4] measured fault zone permeabilities in sandstone outcrops in Arches National Park, Utah, with mini-permeameters. They measure properties normal and parallel to deformation bands and also present a conceptual model of the fault as an anisotropic inhomogeneity embedded in the host rock; they suggest that permeability normal to the fault zone is decreased due to granulation and cataclasis in thin deformation bands that alternate with thin bands of host rock, whereas the permeability parallel to the fault zone is enhanced along the slip plane discontinuity. In other studies, clay smear zones are identified as causing low normal permeabilities while the slip plane discontinuity itself represents a high permeability zone. Bense et al. [9,10] describe the hydraulic properties of faults in unconsolidated sediments in the Roer Valley Rift System in The Netherlands through field and laboratory studies. They conclude that particulate flow in the damage zone surrounding the fault causes the hydraulic conductivity there to be enhanced, while the fault core is characterized by reduced hydraulic conductivity due to clay smearing, grain-scale mixing, and iron-oxide precipitation. Observations of springs along the fault line indicate that vertical groundwater flow along the fault line is enhanced. It appears that the descriptive model of Caine et al. [14] of a fault as an anisotropic inhomogeneity embedded in the host rock is valid for describing the hydrologic properties of faults in a variety of geologic settings. Relating this descriptive model to a quantitative model, however, has proven difficult. It is recognized in the petroleum industry that uncertainty in fault properties is a primary source of uncertainty in reservoir production forecasts [32,24]. Recent studies have focused on quantifying bulk fault properties to be used in regional scale numerical models

1791

of fluid flow [27,40,15,37,19]. Detailed field measurements of fracture geometry [27,19] or stochastic representations of fractures [15,37] are used with numerical simulations to evaluate up-scaled or effective properties of normal fault transmissivity and tangential fault transmissivity. Here, we develop analytical solutions describing groundwater flow through a fault. We model the fault system by its up-scaled properties; we represent the fault as an anisotropic inhomogeneity embedded within the aquifer. Analytical studies of fluid flow through faults are scarce; we have found none that treat the fault in the general form of an anisotropic inhomogeneity. Maas [34] presents the solution for steady flow to a well near a fault in a multiaquifer system; the fault acts to redistribute flow vertically rather than acting as a barrier or conduit for horizontal flow. Shan et al. [39] provide an analytical solution to transient flow near a well in a two-layer aquifer system and propose pumping tests to evaluate fault properties. However, their formulation allows the fault to act only as a source for flow to enter the aquifer, not as a conduit or barrier to horizontal flow; that conceptual model is inconsistent with the anisotropic inhomogeneity described commonly by others. More practical analytical solutions are available in the literature of the petroleum industry: Bixel et al. [11] present the solution for a transient well pumping near a linear inhomogeneity boundary; Yaxley [44] presents a solution for transient flow to a well near an infinite, leaky fault separating two aquifers of different transmissivities. Haneberg [23] considered steady, one-dimensional groundwater flow normal to an idealized fault; he considered the fault to act either as a low or high permeability isotropic inhomogeneity. Using this simple solution, Haneberg [23] and Bense et al. [9] infer fault properties by comparison of head profiles along transects across faults. This approach is useful for estimating the resistance of a fault to normal flow, but provides no estimate of the conductance of the fault to tangential flow. To infer general fault properties, local two-dimensional flow fields near the fault must be investigated. Like Haneburg [23] and Yaxley [44], we consider an idealized fault—vertical, linear, and of infinite extent. Like Yaxley, we include pumping wells in the vicinity and, therefore, non-trivial tangential components of flow are generated along the fault. We limit our analysis to steady flow, and consider the fault transmissivity to be anisotropic with principal directions normal and tangential to the fault. This allows us to model the effects of a general fault as a combined barrier–conduit to flow. 2. Problem description—the fault as an isotropic inhomogeneity The problem domain is illustrated in Fig. 1a; the horizontal, complex z-plane (z = x + iy) is divided into three sub-domains, D1, D*, and D2, with B1 representing the common boundary joining D1 and D*, and B2 joining D* and D2. A pumping well lies at z = zw in D1 and discharges

1792

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

where f(z) is a function that is analytic at z = zw. The boundary conditions along the inhomogeneity boundaries B1 and B2 are given by ( (  U ¼ TT 2 Uþ U ¼ TT 1 Uþ z on B1 z on B ð3Þ 2 W ¼ Wþ W ¼ Wþ

z D2 y

B2 B1

x

D*

where the ‘’ and ‘+’ superscripts indicate evaluation below and above each boundary, respectively. The conditions reflect continuity of head and continuity of the normal component of flow across each boundary. The complex discharge, W(z), is related to the complex potential

D1 zw

(a)

h*

Q

W ¼ T2

T1 x3

T* y

(b) Fig. 1. Definition sketch—Flow to a pumping well near a fault. (a) The horizontal, complex z plane and (b) a section through the fault along the y-axis.

at a steady rate Q [L3/T]; the well is represented by the dot in the figure. Domains D1 and D2 represent the protolith or host aquifer, while D* represents the fault itself. A section through the aquifer along the y-axis is shown in Fig. 1b. The width of the fault zone, domain D*, is h*. We neglect the vertical resistance to flow in the aquifer and adopt the Dupuit approximation [42]. The aquifer transmissivities in domains D1, D*, and D2 are T1, T*, and T2 [L2/T], respectively. We formulate and solve the problem first in terms of a fault with isotropic transmissivity T*; we generalize this solution later to solve the problem of an anisotropic fault with principal transmissivity values of T x and T y . For the case of an isotropic fault, as illustrated in Fig. 1, we seek a solution in the form of a complex potential, X = U + iW, which is an analytic function of the complex coordinate z = x + iy. The discharge potential, U [L3/T], is related to the hydraulic head, / [L], in the aquifer as 8 T /; z in D2 > < 2 ð1Þ U ¼ T  /; z in D > : T 1 /; z in D1 The imaginary part of the complex potential, W [L3/T], is the stream function. The complex potential near the discharge well may be expressed in the form of an expansion of X about z = zw as X¼

Q lnðz  zw Þ þ f ðzÞ 2p

ð2Þ

dX ¼ Qx  iQy dz

ð4Þ

where Qx and Qy [L2/T] are the x and y-components of the discharge vector, respectively. The problem illustrated in Fig. 1 and described above is a degenerate case of the classical problem of flow through two cylindrical inhomogeneities embedded in an infinite domain. The problem is common to several engineering disciplines including groundwater flow, heat conduction, elastostatics, and electrostatics. The solution was originally presented in the field of groundwater mechanics by Kostitsyna [31] as a generalization of Milne-Thomson’s circle theorem [36]. Honein et al. [25] solve the problem of two circular inhomogeneities in terms of anti-plane elastostatics and provide a complete generalization of the circle theorem. In particular, they consider the special case where the two cylinders touch; that solution may be mapped conformally onto the domain in Fig. 1 with a bilinear transformation. Emets and Obnosov [16] and Emets et al. [17] derive and analyze the solution to the problem of two cylindrical inclusions in the context of electrostatics. Kirkham [29] applies complex analysis and presents solutions to similar problems describing groundwater flow to point sinks in layered aquifers, and Maas [35] solves problems in aquifers with an arbitrary number of layers using transform techniques. The solution to the problem illustrated in Fig. 1 may be obtained using the methods described in the references cited above. Alternatively, the solution may be obtained in a convenient form by applying the method of images with respect to the two inhomogeneity boundaries, B1 and B2, as described by Anderson [3]. This approach uses the basic image well solution for a well pumping near a straight inhomogeneity boundary [7, p. 312]. The basic solution is applied iteratively to each boundary, B1 and B2; each iteration results in new image wells that in turn must be imaged about the boundaries. Each additional image well, however, lies outside the problem domain, has a reduced discharge, and lies farther from the real axis than the previous image. The result is a separate expression for the complex potential in each of the three domains. Only the results of this analysis are presented here

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

8 1 n Q P > > > ð1  j1 Þð1 þ j2 Þ 2p ðj1 j2 Þ > > n¼0 > > > > >  lnðz  zw þ 2inh Þ þ TT 21 A; z in D2 > > > > > 1 > > n Q P > ð1  j1 Þ 2p ðj1 j2 Þ lnðz  zw þ 2inh Þ > > > n¼0 > > > > 1 < n Q P ðj1 j2 Þ ð1  j1 Þj2 2p XðzÞ ¼ n¼0 > > >  > >  ln½z  zw  2iðn þ 1Þh  þ TT 1 A; z in D > > > > > > > Q lnðz  zw Þ þ j1 Q lnðz  zw Þ > 2p 2p > > > > 1 P > n > 2 Q > j ð1  j Þ ðj1 j2 Þ 2 > 1 2p > > n¼0 > > :  ln½z  zw  2iðn þ 1Þh  þ A; z in D1

Flow nets featuring contours of hydraulic head (dashed lines) and stream function (solid lines), developed from (5) and (6) are presented in Fig. 3 for two cases. In the first case, the fault has a low transmissivity relative to the adjacent aquifers; in this case, the fault is described commonly as a barrier to flow, or leaky fault. In the second case, the fault has a high transmissivity relative to the adjacent aquifers and the fault is described commonly as a conduit for fluid flow. From the figures we see that both the head and stream function are continuous across both inhomogeneity boundaries. The boundary conditions (3) may be verified visually. In addition, the refraction of streamlines across each inhomogeneity boundary is observed to be correct using the law of refraction of streamlines [7], applied here to boundary B1 as an example tan b T  ¼ tan b1 T 1

ð5Þ where T1  T j1 ¼ T1 þ T

T2  T j2 ¼ T2 þ T

ð7Þ

In (7) b1 is the acute angle that a streamline makes with the y-axis in D1 as it touches the boundary B1, and b* is the corresponding acute angle in domain D*. Further, as the complex potential (5) is an analytic function of the complex variable z, the differential equation is satisfied exactly in each domain.

ð6Þ

and where A is a real constant that may be evaluated from a single reference point of known head in one of the three domains. This solution is depicted graphically in Fig. 2. In the figure, wells are shown as black dots with the strength of the well-displayed next to each dot; only the first three image wells in each infinite series are shown. In the solution valid for z in D1 (Fig. 2a), we see the well of discharge Q (2) in D1 and the first three of an infinite sum of image wells lying outside of D1. In the solution valid for z in D* (Fig. 2b) we see two infinite sums of wells, both lying outside of the domain; in the solution valid for z in D2 (Fig. 2c) we again see an infinite sum of image wells lying outside the domain. In each infinite series, the strength of each image well decreases with distance from the boundaries B1 and B2, and each series converges absolutely [25].

-(1- κ12 ) κ11 κ2 2 Q

1793

3. The fault as an anisotropic inhomogeneity We can generate the solution for the case of an anisotropic inhomogeneity, where the transmissivity of the fault is a tensor with principal axes aligned with the x and yaxes, using standard techniques. The principal values of the fault transmissivity tensor are T x and T y , and the width of the anisotropic fault is w. We transform coordinates in the fault zone, domain D*, following Bear and Dagan [8], such that flow in the transformed domain is governed by the Laplace equation. In D* qffiffiffiffiffiffiffiffiffiffiffiffi ffi the transformation is ~x ¼ x and ~y ¼ y T x =T y . The coordi-

-(1- κ1) κ12 κ2 3 Q 2h* -(1- κ1) κ11 κ2 2 Q

-(1- κ12 ) κ10 κ2 1 Q

zw

zw + κ1 Q

-(1- κ1) κ10 κ2 1 Q y

y x

Q

B1

D*

D2 B2

x

x

B1 (1- κ1) (κ1 κ2) 0 Q

D1

(1- κ1) (1+ κ2 ) ( κ1 κ2 ) 0 Q

zw

zw

y

B2

zw

(1- κ1) (κ1 κ2) 1 Q

(1- κ1) (1+ κ2 ) ( κ1 κ2 ) 1 Q 2h*

(1- κ1) (κ1 κ2) 2 Q

(a)

(b)

(1- κ1) (1+ κ2 ) ( κ1 κ2 ) 2 Q

(c)

Fig. 2. Graphical display of the solution for flow to a well near an isotropic fault. The first three images in the solution valid in each of the three domains (a) D1, (b) D*, and (c) D2 are shown. The dots represent the locations of the well and image wells and the corresponding strength is written adjacent to each well.

1794

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

(a)

(b)

Fig. 3. Example flow fields with the fault represented as an isotropic inhomogeneity for T1/T2 = 2, and jzwj/h* = 2, with (a) T1/T* = 10 and (b) T1/ T* = 0.1. The solid lines are contours of the stream function (DW/Q = 0.05) and the dashed lines are contours of hydraulic head (D(/T1)/Q = 0.05).

nates in D1 remain unaltered, but domain D2 must be shifted upward along the y-axis to account for the transformed width of the domain D*. The transformation of the entire z-plane may be written as hqffiffiffiffiffiffiffiffiffiffiffiffiffi i 8   > z þ iw T  1 ; z in D2 =T > x y > < qffiffiffiffiffiffiffiffiffiffiffiffiffi ð8Þ ~z ¼ ~x þ i~y ¼ x þ iy T  =T  ; z in D > x y > > : z; z in D1 We solve the problem in the transformed domain using the isotropic solution (5), where h* represents the transformed width of domain D*: sffiffiffiffiffi T x  ð9Þ h ¼w T y

8 1 > n Q P > j2 Þ 2p ðf j1 f j2 Þ ð1  f j1 Þð1 þ f > > > n¼0 > > > h qffiffiffiffi i > > > f2 >  ln ez  zew þ 2inw TT x þ TT 21 A; ez in D > y > > > > h qffiffiffiffi i 1 > > n Q P > > ð1  f j1 Þ 2p ðf j1 f j2 Þ ln ez  zew þ 2inw TT x > y > > n¼0 > > > 1 > > n Q P < ð1  f j1 Þf j2 2p ðf j1 f j2 Þ Xð~zÞ ¼ n¼0 > h > qffiffiffiffi i > e > Tx f >  ln e z  z e  2iðn þ 1Þw þ TT 1 A; ez in D > w T y > > > > > Q Q > > lnðez  zew Þ þ f j1 2p lnðez  zew Þ > 2p > > > 1 > P > n Q > f > j2 ð1  f j1 2 Þ 2p ðf j1 f j2 Þ > > > n¼0 > > h qffiffiffiffi i > > > f1 :  ln ez  zew  2iðn þ 1Þw T x þ A; ez in D Ty

The equivalent isotropic transmissivity of the transformed inhomogeneity is, following [8]: qffiffiffiffiffiffiffiffiffiffi Tf ¼ T x T y ð10Þ The coefficients (6) may then be written for the transformed domain, using (10) qffiffiffiffiffiffiffiffiffiffi   T 1  Tf T 1  T x T y qffiffiffiffiffiffiffiffiffiffi f ¼ j1 ¼ ð11Þ T 1 þ Tf T 1 þ T x T y qffiffiffiffiffiffiffiffiffiffi T 2  T x T y qffiffiffiffiffiffiffiffiffiffi f j2 ¼ ð12Þ T 2 þ T x T y The complex potential in the transformed ~z plane is obtained by substituting (9) through (12) into (5), replacing the complex coordinate z with the transformed coordinate ez , and replacing the coefficients j1 and j2 with f j1 and f j2 . We obtain

ð13Þ f

f1 , D , and D f2 are the transformed domains in the ~z where D plane. Substitution of (8) for ~z in (13) yields the solution in the physical, z-plane, where we note that in each expression (13) zew ¼ zw and zew ¼ zw . In the physical plane, a complex potential exists only in domains D1 and D2 where the flow is irrotational. We obtain 8 1 n > Q P > f f ð1  j Þð1 þ j Þ ðf j1 f j2 Þ ln ½z  zw > 1 2 2p > > n¼0 > > i qffiffiffiffi > > Tx > >  iw þ TT 21 A; z in D2 þ ið2n þ 1Þw  > Ty > < Q Q XðzÞ ¼ 2p lnðz  zw Þ þ f j1 2p lnðz  zw Þ > > 1 > P > Q > j1 2 Þ 2p ðf j1 f j 2 Þn f j2 ð1  f > > > n¼0 > > h qffiffiffiffi i > > > :  ln z  zw  2iðn þ 1Þw T x þ A; z in D1 Ty

ð14Þ

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

(a)

1795

(b)

Fig. 4. Example flow fields for an anisotropic fault. In both figures T1/T2 = 2, jzwj/w = 2 with (a) T 1 =T x ¼ 10, T x =T y ¼ 0:01 and (b) T 1 =T x ¼ 0:1, T x =T y ¼ 100. The solid lines are contours of the stream function (DW/Q = 0.05) and the dashed lines are contours of hydraulic head (D(/T1)/Q = 0.05).

In domain D*, where the flow is rotational, we write two real functions—one for evaluating the hydraulic head and one for evaluating the stream function: ( /ðx; yÞ ¼ RX½ez ðx; yÞ= Tf  z in D ð15Þ Wðx; yÞ ¼ IX½ez ðx; yÞ where ~zðx; yÞ is given by (8). Eqs. (14) and (15) provide the solution to the problem when the fault is represented as an anisotropic inhomogeneity of width w. Flow nets for two cases where the fault is represented as an anisotropic inhomogeneity are provided in Fig. 4. In Fig. 4a, the anisotropy ratio of the fault, T x =T y , is 0.01 and in Fig. 4b, T x =T y ¼ 100. Interpretation of flow based on head patterns is difficult in the anisotropic domain as streamlines are not normal to head contours. For example, in Fig. 4a the head gradient within the fault is directed primarily tangential to the fault, and there is little change in head across the fault; based on head contours alone, this seems to imply that the fault behaves as a conduit with a significant tangential component of flow within the fault. The streamlines, however, indicate that this is not the case; the anisotropy limits the tangential component of flow within the fault and forces the flow to be primarily normal to the fault. The streamlines within the fault and in D1 indicate the behavior is that of a leaky fault rather than a conductive fault. The opposite effect is seen in Fig. 4b where the head gradient within the fault implies that flow is primarily normal to the fault, while the streamlines show that flow within the fault is primarily tangential to the fault. 4. Fault conductance and resistance In regional models of groundwater flow, a fault is often viewed as a linear feature whose width is negligible when compared to the scale of the model. We will show here that the effects of an anisotropic fault on flow in the adjacent aquifers is independent of the width of the fault; the prop-

erties of the fault may be described exactly by two lumped parameters—fault resistance and fault conductance. We may eliminate domain D*, which is thin with respect to the scale of the model, while including the effects of the fault through the lumped parameters. We define the fault conductance, Cf [L3/T], to be the product of the fault width w and the fault transmissivity tangential to the fault, T x C f ¼ wT x

ð16Þ

In a similar fashion, we define the fault resistance, Rf [T/L], to be the fault width divided by the fault transmissivity normal to the fault, T y w Rf ¼  ð17Þ Ty We consider the limiting case of flow in the fault where the thickness of the fault w and transmissivity in the normal direction T y vanish, while holding their ratio constant. At the same time we let the fault transmissivity in the tangential direction T x become infinite while the product of w and T x remains constant. Thus, we apply the limit to the complex potential (14) for w ! 0, T y ! 0, and T x ! 1 while holding w=T y ¼ Rf and wT x ¼ C f . To simplify notation, we express this limit as limRf ;Cf . We note the following: sffiffiffiffiffi T x pffiffiffiffiffiffiffiffiffiffi ¼ C f Rf lim w ð18Þ Rf ;C f T y and, from (12) and (11) pffiffiffiffiffiffiffiffiffiffiffiffi T 1  C f =Rf pffiffiffiffiffiffiffiffiffiffiffiffi ¼ K 1 f lim j1 ¼ Rf ;C f T 1 þ C f =Rf pffiffiffiffiffiffiffiffiffiffiffiffi T 2  C f =Rf pffiffiffiffiffiffiffiffiffiffiffiffi ¼ K 2 j2 ¼ lim f Rf ;C f T 2 þ C f =Rf

ð19Þ ð20Þ

We apply the limit to the complex potential in D1 and D2, (14), using (18)–(20) to obtain

1796

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

-(1-K12 ) K11K2 2 Q 2 Cf Rf

-(1-K12 )

K10K2 1 Q

D2

zw + K1 Q

y

y

x

x

Q zw

D1

zw -i Cf Rf +(1-K1) (1+K2 ) ( K1K2 ) 0 Q +(1-K1) (1+K2 ) ( K1K2 ) 1 Q 2 Cf Rf +(1-K1) (1+K2 ) ( K1K2 ) 2 Q

Fig. 5. Graphical display of the solution for flow to a well near a general fault boundary.

(a)

(b)

Fig. 6. Example flow fields for a general fault characterized by a resistance and a conductance. In both figures T1/T2 = 2 with (a) RfT1/jzwj = Cf/ (T1jzwj) = 0.01 and (b), RfT1/jzwj = Cf/(T1jzwj) = 10. The solid lines are contours of the stream function (DW/Q = 0.05) and the dashed lines are contours of hydraulic head (D(/T1)/Q = 0.05).

8 1 n Q P > ðK 1 K 2 Þ ln½z  zw ð1  K 1 Þð1 þ K 2 Þ 2p > > > n¼0 > > pffiffiffiffiffiffiffiffiffiffi > > > þ ið2n þ 1Þ C f Rf Þ þ TT 2 A; z in D2 > 1 < Q Q lim XðzÞ ¼ 2p lnðz  z Þ þ K lnðz  zw Þ w 1 2p Rf ;C f > > > 1 P > n 2 Q > > > >  K 2 ð1  K 1 Þ 2p n¼0ðK 1 K 2 Þ ln½z  zw > > : pffiffiffiffiffiffiffiffiffiffi  2iðn þ 1Þ C f Rf Þ þ A; z in D1

tance, using (17) and (16). The flow fields in D1 and D2 presented in Figs. 4a and 6a are identical, as are the flow fields in D1 and D2 of Figs. 4b and 6b, providing visual evidence that the flow in D1 and D2 is independent of the fault width, w. We note in Fig. 6 that both the head and stream function are discontinuous across the fault, although the jump in head in Fig. 6a is small. ð21Þ

The solution (21) is displayed graphically in Fig. 5. When taking the limit, the domain D* vanishes, although thep trans-ffi f remains of constant width h ¼ ffiffiffiffiffiffiffiffiffi formed domain D C f Rf f (13) remains un(see (9) and (18)), and the flow field in D changed. The effects of the fault on flow in the adjacent aquifer is dependent only on the two fault parameters, Rf and Cf, and is independent of the physical width, w, of the fault. Example flow nets generated from (21) are given in Fig. 6. The two cases presented in Fig. 6 correspond to the cases presented in Fig. 4, where the fault transmissivities and width are replaced by fault resistance and conduc-

5. Limiting cases The solution (21) consists pffiffiffiffiffiffiffiffiffiffi of infinite sums of wells spaced at intervals of 2 C f Rf and lying along the imaginary axis, as illustrated in Fig. 5. When either parameter, Rf or Cf, is small, the distance between image wells is small. We consider the special cases where the fault acts only as a leaky barrier to normal flow (Cf = 0), and where the fault acts only as a conduit to tangential flow (Rf = 0). These two limiting cases are the end-member fault types, identified by Sigda et al. [41], of the qualitative model of fault properties presented by Caine et al. [14]. In both cases the distance between image wells vanishes and the infinite

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

sums of wells in the solution (21) become continuous distributions of wells along semi-infinite lines. For convenience in evaluating these special cases we express the infinite sums in (21) as discrete line sinks, following [3]. For example, for z in D2 we write XðzÞ ¼ ð1  K 1 Þð1 þ K 2 Þ

Q 2p

n¼0 1 X T2 rm T2 lnðz  dm Þ Dn þ A A¼ T1 2p T1 m¼1

ð22Þ

where m = n + 1 and dm is the location of well m, defined as pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi dm ¼ zw  ið2n þ 1Þ C f Rf ¼ zw  ið2m  1Þ C f Rf ð23Þ The discharge of well m is rm Dn where Dn represents the distance pffiffiffiffiffiffiffiffiffiffi between well m and well m + 1, and is equal to 2 C f Rf . From (22) Q ð1  K 1 Þð1 þ K 2 Þ m rm ¼ pffiffiffiffiffiffiffiffiffiffi ðK 1 K 2 Þ K 1K 2 2 C f Rf

ð24Þ

We express the integer m in terms of the location of well m by inverting (23) pffiffiffiffiffiffiffiffiffiffi z w þ i C f Rf  dm pffiffiffiffiffiffiffiffiffiffi m¼ ð25Þ 2i C f Rf By substitution of (25) for m into (24) we obtain the strength of well m as a function of its location rm Dn ¼

¼

Qð1  K 1 Þð1 þ K 2 Þ pffiffiffiffiffiffiffiffiffiffi 2 C f Rf K 1 K 2 pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi  ðK 1 K 2 Þðzw þi Cf Rf dm Þ=½2i Cf Rf  Dn Qð1  K 1 Þð1 þ K 2 Þ pffiffiffiffiffiffiffiffiffiffi 2 C f Rf K 1 K 2   pffiffiffiffiffiffiffiffiffiffi lnðK 1 K 2 Þ pffiffiffiffiffiffiffiffiffiffi ðzw þ i C f Rf  dm Þ Dn  exp 2i C f Rf

ð26Þ

Finally, combination of (22) and (26) gives the complex potential for z in D2 X¼

Qð1  K 1 Þð1 þ K 2 Þ pffiffiffiffiffiffiffiffiffiffi 4p C f Rf K 1 K 2   1 X pffiffiffiffiffiffiffiffiffiffi lnðK 1 K 2 Þ T2 pffiffiffiffiffiffiffiffiffiffi ðzw þ i C f Rf  dm Þ lnðz  dm Þ Dn þ A  exp T1 2i C f Rf m¼1 ð27Þ

where dm is given by (23). In a similar fashion, for z in D1, the complex potential (21) may also be written using the discrete form of a line sink Q Q Qð1  K 2 Þ lnðz  zw Þ þ K 1 lnðz  zw Þ  pffiffiffiffiffiffiffiffiffiffi1 2p 2p 4p C f Rf K 1   1 X lnðK 1 K 2 Þ pffiffiffiffiffiffiffiffiffiffi ðdm  zw Þ lnðz  dm ÞDn þ A  exp 2i C f Rf m¼1

where pffiffiffiffiffiffiffiffiffiffi dm ¼ zw þ 2im C f Rf

ð29Þ

Eqs. (27) and (28) will be used to evaluate the complex potential for the two limiting cases, Cf = 0 and Rf = 0. 5.1. The leaky fault

1 X pffiffiffiffiffiffiffiffiffiffi  ðK 1 K 2 Þn ln½z  zw þ ið2n þ 1Þ C f Rf Þ

þ

1797

X¼

ð28Þ

When the conductance of the fault to tangential flow vanishes, the result is a purely leaky fault. This leaky fault is a limiting case of the anisotropic fault that allows for a one-dimensional shear flow (Qy(x)) to develop: flow is one-dimensional and normal to the fault (in the y-direction) everywhere in the fault domain, but varies in magnitude along the fault. As flow in the domain D* is one-dimensional, the domain may be replaced exactly by internal boundary conditions. For a leaky fault, the appropriate conditions joining domains D1 and D2 along the real axis are [42] W ¼ Wþ þ Q y ¼ Qy ¼





þ

1 U U  Rf T 1 T2

ð30Þ ð31Þ

The first condition (30) reflects a continuous normal component of flow across the fault. The second condition (31) states that the normal component of flow is equal to the jump in head across the fault divided by the fault resistance. We generate the solution to the problem of a leaky fault by considering it to be a limiting case of an anisotropic inhomogeneity. We consider the coefficients appearing in the complex potential, (27) and (28). The first coefficient in (27) may be expanded using (19) and (20) as   Qð1  K 1 Þð1 þ K 2 Þ Q T1 pffiffiffiffiffiffiffiffiffiffi ð32Þ ¼ 2 pRf T 1  ðC f =Rf Þ 4p C f Rf K 1 K 2 Taking the limit for Cf ! 0 yields lim

C f !0

Qð1  K 1 Þð1 þ K 2 Þ Q pffiffiffiffiffiffiffiffiffiffi ¼ pT 1 Rf 4p C f Rf K 1 K 2

ð33Þ

In a similar fashion the limit of the coefficient in (28) may be obtained lim

C f !0

Qð1  K 21 Þ Q pffiffiffiffiffiffiffiffiffiffi ¼ 4p C f Rf K 1 pT 1 Rf

ð34Þ

Finally, in both expressions (27) and (28), the limit of the coefficient within the argument of the exponential must be evaluated. We expand the logarithmic term in a Taylor series about Cf = 0 to obtain   lnðK 1 K 2 Þ i 1 1 lim pffiffiffiffiffiffiffiffiffiffi ¼ þ ð35Þ C f !0 2i C f Rf Rf T 1 T 2 In the limit for Cf ! 0, the discrete sum of point sinks in (27) becomes a continuous distribution of sinks along the line extending from zw to zw  i1. The complex potential for z in D2 becomes, using (33) and (35)

1798

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

   i 1 1 exp  þ ðd  zw Þ Rf T 1 T 2 zw T2  lnðz  dÞ dd þ A; z in D2 ð36Þ T1

Q i lim X ¼ C f !0 p T 1 Rf

Z

zw i1



where the complex number d is the integration variable that denotes a point of the line extending along the imaginary axis from zw to negative infinity and may be written in terms of n as d ¼ zw  in;

06n61

We make the following change of variables:   i 1 1 þ Z2 ¼  ðz  zw Þ Rf T 1 T 2   i 1 1 þ D¼ ðd  zw Þ Rf T 1 T 2 and the complex potential becomes  Z 1 Q T2 X¼ expðDÞ lnðZ 2  DÞ dD p T1 þ T2 0     Z 1 i 1 1 T2  ln  þ expðDÞ dD þ A Rf T 1 T 2 T1 0

ð37Þ

ð38Þ ð39Þ

ð40Þ

The first integral may be integrated by parts and the second may be integrated directly to obtain   Q T2 X¼ lnðz  zw Þ p T1 þ T2   Q T2 T2 þ expðZ 2 ÞE1 ðZ 2 Þ þ A; z in D2 ð41Þ p T1 þ T2 T1 where E1(z) is the exponential integral [1]. Some properties of the exponential integral and the function exp(z)E1(z) are provided in Appendix A, and Appendix B provides expressions for the complex discharge W for both of the limiting solutions. A similar procedure may be taken to evaluate the complex potential for z in D1 when we take the limit of (28) for Cf ! 0. We obtain the following expression for z in D1, using (34) and (35):   Q Q T1  T2 lnðz  zw Þ þ X¼ lnðz  zw Þ 2p 2p T 1 þ T 2   Q T2  ð42Þ expðZ 1 ÞE1 ðZ 1 Þ þ A; z in D1 p T1 þ T2 where Z1 ¼

i Rf



 1 1 þ ðz  zw Þ T1 T2

ð43Þ

We verify in Appendix C that the complex potential, (41) and (42), satisfies the boundary conditions (30) and (31) exactly. The solution, (41) and (42), is the steady-state equivalent of the solution to the transient problem solved by Yaxley [44]. 5.2. The conductive fault When the resistance of the fault to normal flow vanishes, the result is a conductive fault. This also is a limiting type

of anisotropy, equivalent to Dupuit–Forchheimer flow, where the head is one-dimensional in the fault (/(x)), although the stream function is two-dimensional. In this case, the domain D* may be replaced exactly by the following boundary conditions joining domains D1 and D2 along the real axis: U Uþ ¼ ð44Þ T1 T2 Q Qþ C f x ¼ C f x ¼ W  Wþ ð45Þ T1 T2 The first condition (44) ensures that the head is continuous across the fault; the second (45) states that the flow within the fault is equal to the jump in stream function across the fault. The complex potential that satisfies (44) and (45) may be found by following a limiting process similar to that taken for the leaky fault. We expand the coefficients that appear in the complex potential, (27) and (28) and take the limit for Rf ! 0. Only the results of this analysis are presented here, although we verify that the solution satisfies the boundary conditions, (44) and (45) in Appendix D. The complex potential for the conductive fault is 8   Q T2 > > > lnðz  zw Þ > > p T1 þ T2 > >   > > > Q T2 T2 > > < þ p T þ T expðZ 2 ÞE1 ðZ 2 Þ þ T A; z in D2 1 2 1   X¼ > Q Q T  T 1 2 > > lnðz  zw Þ þ lnðz  zw Þ > > 2p 2p T 1 þ T 2 > >   > > > Q T1 > > z in D1 expðZ 1 ÞE1 ðZ 1 Þ þ A; : þ p T1 þ T2 ð46Þ where i ðT 1 þ T 2 Þðz  zw Þ Cf i Z 2 ¼  ðT 1 þ T 2 Þðz  zw Þ Cf

Z1 ¼

ð47Þ

6. Discussion: the conductive fault as a barrier to flow Examples of flow nets for the two limiting cases and for a general fault of zero width are presented in Fig. 7. The three fault types are distinguished by the behavior of the head or stream function across the fault. A leaky fault is depicted in Fig. 7a; the stream function is continuous while the head jumps across the fault. A conductive fault is depicted in Fig. 7b; here the head is continuous across the fault, but the stream function jumps. Finally, a general fault is depicted in Fig. 7c; in the case of the general fault, both the head and stream function jump across the fault. It is of interest to note that the form of the complex potentials for the two limiting cases are similar to a solution presented by Bakker and Anderson [6] to an unrelated problem with internal boundaries. In that problem, the internal boundary represents a stream removing water from the aquifer

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

(a)

1799

(b)

(c) Fig. 7. Example flow fields for T1/T2 = 1. A leaky fault is shown in (a) with RfT2/jzwj = 5. A conductive fault is shown in (b) with Cf/(T1jzwj) = 5. A general fault is shown in (c) with RfT2/jzwj = Cf/(T1jzwj) = 5. The solid lines are contours of the stream function (DW/Q = 0.05) and the dashed lines are contours of hydraulic head (D(/T1)/Q = 0.05).

and reflects a jump in the imaginary part of the derivative of the complex potential rather than a jump in the complex potential itself. A comparison of the complex potentials for the leaky fault (41) and the conductive fault (46) shows that the expressions are similar in form. In particular, the expressions valid for z in D2 are identical except for the form of the variable Z2. For the leaky fault we have (38) and for the conductive fault (47). The aquifer transmissivities for the leaky fault are combined in parallel in the expression for Z2, while the transmissivities for the conductive fault are combined in series. Setting these two expressions equal, we find that the flow fields in D2 are identical for a conductive fault and a leaky fault when the following condition is satisfied: Cf ¼ Rf T 2 T1

ð48Þ

This behavior may be observed in Fig. 7a and b for the case that T1/T2 = 1; in Fig. 7a, RfT2/jzwj = 5 and in Fig. 7b, Cf/ (T1jzwj) = 5, which satisfy condition (48). The flow fields in D2 are identical although their behavior in D1 is quite different.

This same behavior can be observed in the potential where the fault is modeled as an isotropic inhomogeneity (5) when comparing results for a high-transmissivity fault with results for a low transmissivity fault. We find for T1/T2 = 1 that the flow fields are identical in D2 for T*/ T1 = a (where a is a constant) and T*/T1 = 1/a. This may be seen by writing the strength of the nth image well for z in D2, Qn, in terms of the dimensionless parameters T1/ T* and T1/T2. From (5), the strength of the nth image well is given by Qn ¼ Qð1  j1 Þð1 þ j2 Þðj1 j2 Þ

n

ð49Þ

where Q is the strength of the pumping well in D1. We use (6) to write this in terms of the dimensionless parameters T1/T* and T1/T2 to obtain h  in    T 1  1 T 1 T 2  1 T T T1 T1 T2 ð50Þ Qn ¼ 4Q  h inþ1

T T T T 1 T 1 1 2 þ 1 þ 1   T T T1 We write this in an alternate form by multiplying both the numerator and denominator in (50) by (T*/T1)2(n+1) to obtain

1800

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

h  in       1  TT 1 TT 21  TT 1 T T2 Qn ¼ 4Q h  inþ1   T1 T1 1 þ TT 1 TT 21 þ TT 1

ð51Þ

For the special case T1/T2 = 1, we obtain from (50) and (51), respectively  2n

2n   T 1   T   1 T1 T1 T   1 Qn ¼ 4Q  T 2ðnþ1Þ ¼ 4Q  2ðnþ1Þ T  1 T T1 T þ1 þ 1 T T1 ð52Þ The transmissivity T*/T1 may be replaced by T1/T* in the expression for Qn without affecting the value of Qn; for example, the complex potentials (5) for T*/T1 = a and T*/T1 = 1/a are equivalent in D2. This analysis shows that inferring the properties of a fault based on measured flow or heads in the domain opposite the fault from the feature causing flow—in this case the well—is not possible. This also has implications for the effectiveness of engineered barriers such as slurry walls and high permeability walls: both may act as barriers to flow. In the setting of a faulted aquifer, the conductive fault is as effective as the leaky fault at blocking groundwater in D2 from entering D1. However, the path taken, the head loss, and thus the travel time through the fault are significantly different. Fig. 8 shows profiles of the potential along the y-axis for the three flow fields presented in Fig. 7. As discussed earlier, the profile for the leaky and conductive faults are identical in D2. We also see that the profile for the general fault is very similar in D2 as well as being similar to the profile for the leaky fault in D1. The jump in the dimensionless discharge potential across the fault for the general fault is only about 14% less than the case of the leaky fault, which has the same resistance; although the profile of the discharge potential for the two cases are similar both in D1 and D2, the flow fields illustrated in Fig. 7a and c are significantly different due to the high conductance of the fault in the general case. While it may be possible to use the observed head change across a fault to estimate the resistance of the

0.4 0.2

Φ/Q

0.0

General Fault

-0.2 Conductive Fault

-0.4

7. Conclusions Explicit analytic expressions describing groundwater flow to a well through a fault have been presented for several common descriptions of the fault zone. The most general solution represents the fault zone as an anisotropic inhomogeneity, allowing the fault to act as a combined barrier–conduit to fluid flow, as is commonly described in the literature. The effects of an anisotropic fault on flow in the adjacent aquifers can be characterized by two parameters—fault resistance, Rf, and fault conductance, Cf—and are independent of the width of the fault. When one of these parameters is zero, either the discharge or the head distribution within the fault is one-dimensional, and the effects of the fault may be incorporated exactly by replacing the fault domain with internal boundary conditions. In the general case, where both parameters are non-zero, internal boundary conditions cannot be used to replace the anisotropic domain. Although the physical fault width w is allowed to vanish, the width of the transformed domain h* remains unchanged. Both discharges and head distributions are two-dimensional within the fault and the domain D* cannot be replaced exactly by internal boundary conditions. Using the analytic solutions, we have found that the term ‘‘barrier’’, commonly used to describe a leaky fault, is inappropriate. We have shown that a leaky fault and a conductive fault both act as barriers to flow in the domain opposite the pumping well; the flow fields for both types of faults were shown to be identical in many cases. This result is not intuitive, although similar behavior has been observed in the context of engineered barriers for waste isolation [28]. The analytic solutions developed here can be used, for example, to design steady-state pumping tests to infer fault properties based on observations of drawdown and well discharge. Field observations of head have been used commonly to infer fault properties, both qualitatively and quantitatively [30,23,10]. We have shown, however, that interpretation of fault properties is difficult, even under the controlled conditions represented by the analytical solutions. The effect of fault resistance is easily observed as a jump in head across the fault; fault conductance, however, is difficult to infer from observed head patterns.

Leaky Fault

-0.6 -0.8 -2.0

fault, it would be difficult to distinguish between cases of a fault with no conductance (the leaky fault) and a general fault with high conductance, even in cases when a well is used to produce a significant tangential gradient along the fault.

Acknowledgments -1.0

0.0

1.0

2.0

y / |zw|

Fig. 8. A profile of the potential along the y-axis comparing the behavior of the leaky fault, the conductive fault and the general fault. The thick line represents the general fault.

Portions of this work were supported by a grant from the Environmental Research Initiative Committee of the University of South Carolina. Elizabeth Mesa assisted with preparation of the figures.

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

1801

Appendix A. Properties of exp(z)E1(z)

Appendix C. Boundary conditions: the leaky fault

The exponential integral is defined by Abromowitz and Stegun [1] as Z 1 expðtÞ dt ðA:1Þ E1 ðzÞ ¼ t z

For the leaky fault, the internal boundary conditions are given by (30) and (31). We note that

and the first derivative is dE1 ðzÞ expðzÞ ¼ dz z

ðA:2Þ

The function exp(z)E1(z) represents a semi-infinite line dipole of exponentially decaying strength distribution, as discussed by Anderson [3]. The derivative of this function is d 1 ½expðzÞE1 ðzÞ ¼  þ expðzÞE1 ðzÞ dz z

Iflnðx  zw Þg ¼ Iflnðx  zw Þg ¼ argðx  zw Þ

Taking the imaginary part of the complex potential (42) for z in D1, and evaluating along y = 0 using (C.1) we obtain   Q T2 W ¼ argðx  zw Þ p T1 þ T2   Q T2   ðC:2Þ IfexpðZ  1 ÞE 1 ðZ 1 Þg p T1 þ T2 where

ðA:3Þ

Both the exponential function and the exponential integral satisfy symmetry relationships [1]

Z 1

i ¼ Rf



 1 1 þ ðx  zw Þ T1 T2

From (A.4) and (A.5) we can evaluate the symmetry properties of the semi-infinite line dipole:

For z in D2 (41) we obtain   Q T2 þ W ¼ argðx  zw Þ p T1 þ T2   Q T2 þ þ IfexpðZ þ 2 ÞE 1 ðZ 2 Þg p T1 þ T2

expðzÞE1 ðzÞ ¼ expðzÞE1 ðzÞ

where

expðzÞ ¼ expðzÞ

ðA:4Þ

E1 ðzÞ ¼ E1 ðzÞ

ðA:5Þ

ðA:6Þ

i ¼ Rf

 1 1 þ ðx  zw Þ T1 T2

Zþ 2

RfexpðzÞE1 ðzÞg ¼ RfexpðzÞE1 ðzÞg

ðA:7Þ

We note from (C.3) and (C.5) that

IfexpðzÞE1 ðzÞg ¼ IfexpðzÞE1 ðzÞg

ðA:8Þ

Appendix B. The complex discharge We provide explicit expressions for the complex discharge, W (4), for the two limiting cases presented. For the leaky fault, we find using (4), (A.3), (41)–(43) 8 iQ > > expðZ 2 ÞE1 ðZ 2 Þ; z in D2 > > > pRf T 1 > >   < Q 1 1 ðB:1Þ W ðzÞ ¼  þ > 2p z  z z  zw w > > > > iQ > > : þ expðZ 1 ÞE1 ðZ 1 Þ; z in D1 pT 1 Rf where Z2 and Z1 are given by (38) and (43), respectively. For the conductive fault we obtain, using (4), (A.3), (46) and (47) 8 iQT 2 > > expðZ 2 ÞE1 ðZ 2 Þ; z in D2 > > > pC f > >     < Q 1 T2 1 W ðzÞ ¼   > 2p z  z z  T þ T zw w 1 2 > > > > iQT 1 > > :  expðZ 1 ÞE1 ðZ 1 Þ; z in D1 pC f ðB:2Þ

ðC:3Þ

ðC:4Þ



Separating the left- and right-hand sides of (A.6) into real and imaginary parts yields the following relationships:

Here, Z1 and Z2 are given by (47).

ðC:1Þ

þ Z 1 ¼ Z2

ðC:5Þ

ðC:6Þ

and from the symmetry relationships (A.8), we see that the boundary condition (30) is satisfied; the normal component of flow is continuous across the fault, as may be seen by comparing (C.2) and (C.4). Next we check the jump in head across the fault. Using the relationship Rflnðx  zw Þg ¼ Rflnðx  zw Þg ¼ ln jx  zw j and using (41) and (42), we find for z in D1   U Q 1 ¼ ln jx  zw j p T1 þ T2 T1   Q T2 A   RfexpðZ  1 ÞE 1 ðZ 1 Þg þ pT 1 T 1 þ T 2 T1 Similarly, we find for z in D2   Uþ Q 1 ¼ ln jx  zw j p T1 þ T2 T2   Q 1 A þ þ RfexpðZ þ 2 ÞE 1 ðZ 2 Þg þ p T1 þ T2 T1 Using the symmetry relationship (A.7) we obtain   1 U Uþ Q þ  RfexpðZ þ ¼ 2 ÞE 1 ðZ 2 Þg Rf T 1 pT 1 Rf T2

ðC:7Þ

ðC:8Þ

ðC:9Þ

ðC:10Þ

1802

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

From (B.1) we evaluate Qþ y as Q þ þ Qþ RfexpðZ þ y ¼ IðW Þ ¼  2 ÞE1 ðZ 2 Þg pT 1 Rf

References ðC:11Þ

By comparison of (C.10) and (C.11) we see that the second boundary condition (31) is also satisfied. Appendix D. Boundary conditions: the conductive fault For the conductive fault, the internal boundary conditions are given by (44) and (45). We take the real part of the complex potential (46) valid for z in D1, and evaluate it along the real axis to obtain   U Q 1 ¼ ln jx  zw j p T1 þ T2 T1   Q 1 A  ðD:1Þ þ RfexpðZ  1 ÞE 1 ðZ 1 Þg þ p T1 þ T2 T1 From the expression valid for z in D2 we obtain   Uþ Q 1 ¼ ln jx  zw j p T1 þ T2 T2   Q 1 A þ þ RfexpðZ þ 2 ÞE 1 ðZ 2 Þg þ p T1 þ T2 T1

ðD:2Þ

where use is made of (C.7). Using (A.7) and (A.8), we see that the expressions (D.1) and (D.2) are equivalent and the boundary condition (44) is satisfied. To check the second boundary condition (45), we evaluate the imaginary part of the complex potential along the real axis. From the expression valid for z in D1 we obtain   Q T2  W ¼ argðx  zw Þ p T1 þ T2   Q T1  ðD:3Þ þ IfexpðZ  1 ÞE1 ðZ 1 Þg p T1 þ T2 where use is made of (C.1). Similarly, from the expression valid for z in D2 we obtain   Q T2 þ W ¼ argðx  zw Þ p T1 þ T2   Q T2   ðD:4Þ IfexpðZ  1 ÞE 1 ðZ 1 Þg p T1 þ T2 where use is made of (A.8) and (C.6). Subtracting (D.4) from (D.3), we find Q  W  Wþ ¼ IfexpðZ  1 ÞE 1 ðZ 1 Þg p Q þ ¼  IfexpðZ þ ðD:5Þ 2 ÞE 1 ðZ 2 Þg p From (B.2) for z = x in D2 we have Cf þ Q þ Qx ¼  IfexpðZ þ ðD:6Þ 2 ÞE1 ðZ 2 Þg p T2 By comparison of the two expressions (D.5) and (D.6) we find that the second boundary condition (45) is also satisfied exactly.

[1] Abramowitz M, Stegun IA. Handbook of mathematical functions with formulas, graphs, and mathematical tables. New York, NY: Dover; 1965. [2] Ahlbom K, Smellie JAT. Overview of fracture zone project at Finnsjon, Sweden. J Hydrol 1991;126:1–15. [3] Anderson EI. The method of images for leaky boundaries. Adv Water Resour 2000;23:461–74. [4] Antonellini M, Aydin A. Effect of faulting on fluid flow in porous sandstones: petrophysical properties. AAPG Bull 1994;78:355– 77. [5] Babiker M, Gudmundsson A. The effects of dykes and faults on groundwater flow in an arid land: the Red Sea Hills, Sudan. J Hydrol 2004;297:256–73. [6] Bakker M, Anderson EI. Steady flow to a well near stream with a leaky bed. Ground Water 2003;41(6):833–40. [7] Bear J. Dynamics of fluids in porous media. New York, NY: Dover Publications, Inc; 1988. [8] Bear J, Dagan G. The relationship between solutions of flow problems in isotropic and anisotropic soils. J Hydrol 1965;3:88–96. [9] Bense VF, van Balen RT, De Vries JJ. The impact of faults on the hydrogeological conditions in the Roer Valley Rift System: an overview. Netherlands J Geosci 2003;82(1):317–20. [10] Bense VF, Van den Berg EH, Van Balen RT. Deformation mechanisms and hydraulic properties of fault zones in unconsolidated sediments: the Roer Valley Rift System, the Netherlands. Hydrogeol J 2003;11(3):319–32. [11] Bixel HC, Larkin BK, Van Poolen HK. Effect of linear discontinuities on pressure build-up and drawdown behavior. J Petrol Technol 1963;15:885–95. [12] Bredehoeft JD, Belitz K, Sharp-Hansen S. The hydrodynamics of the Big Horn Basin: a study of the role of faults. AAPG Bull 1992;76(1):530–46. [13] Bredehoeft JD. Fault permeability near Yucca Mountain. Water Resour Res 1997;33(11):2459–63. [14] Caine JC, Evans JP, Forster CB. Fault zone architecture and permeability structure. Geology 1996;24(11):1025–8. [15] Caine JS, Forester CB. Fault zone architecture and fluid flow: insights from field data and numerical modeling. In: Haneberg WC et al., editors. Faults and subsurface fluid flow in the shallow crust. Geophysical monograph series, vol. 113. Washington, DC: AGU; 1999. p. 101–27. [16] Emets YuP, Obnosov YuV. An accurately solvable problem of the mutual effect of inclusions in the theory of heterogeneous media. J Appl Mech Tech Phys 1990;31:21–9. [17] Emets YuP, Obnosov YuV, Onofriichuk YuP. Interaction between touching circular dielectric cylinders in a uniform electric field. J Tech Phys 1993;38(12):1043–7. [18] Evans JP, Forster CB, Goddard JV. Permeability of fault-related rocks, and implications for hydraulic structure of fault zones. J Struct Geol 1997;19:1393–404. [19] Flodin EA, Durlofsky LJ, Aydin A. Upscaled models of flow and transport in faulted sandstone: boundary condition effects and explicit fracture modelling. Petrol Geosci 2004;10:173–81. [20] Forster CB, Evans JP. Hydrogeology of thrust faults and crystalline thrust sheets: results of combined field and modeling studies. Geophys Res Lett 1991;18(5):979–82. [21] Fowles J, Burley S. Textural and permeability characteristics of faulted, high porosity sandstones. Marine Petrol Geol 1994;11: 608–23. [22] Ganser DR. Hydrogeologic characteristics of the Ramapo Fault, Northern New Jersey. Ground Water 1988;25(6):664–71. [23] Haneberg WC. Steady state groundwater flow across idealized faults. Water Resour Res 1995;31(7):1815–20. [24] Harris D, Yielding G, Levine P, Maxwell G, Rose PT, Nell P. Using shale gouge ratio (SGR) to model faults as transmissibility barriers in

E.I. Anderson / Advances in Water Resources 29 (2006) 1790–1803

[25] [26]

[27]

[28]

[29] [30]

[31]

[32]

[33]

[34]

reservoirs: an example from the Strathspey Field, North Sea. Petrol Geosci 2002;8:167–76. Honein E, Honein T, Herrmann G. On two circular inclusions in harmonic problems. Quart Appl Math 1992;1(3):479–99. Huntoon PW, Lundy DA. Fracture-controlled ground-water circulation and well sitting in the vicinity of Laramie, Wyoming. Ground Water 1979;17(5):463–9. Jourde H, Flodin EA, Aydin A, Durlofsky LJ, Wen XH. Computing permeability of fault zones in eolian sandstone from outcrop measurements. AAPG Bull 2002;86(7):1187–200. Kacimov AR, Obnosov YuV. Minimization of ground water contamination by lining of a porous waste repository. Proc Indian Natl Sci Acad, Part A 1994;60(6):783–92. Kirkham D. Seepage of artesian and surface water into drain tubes in stratified soil. Trans AGU 1954;35(5):775–90. Kolm KE, Downey JS. Diverse flow patterns in the aquifers of the Amargosa Desert and vicinity, Southern Nevada and California. Bull Assoc Eng Geol 1994;31(1):33–47. Kostitsyna LI. Translational filtration flow across a system of spherical layers with different permeabilities. In: Hydromechanics, Collection of Papers, 1973 [in Russian]. Publishing house of the NK Krupskaya Moscow Region Teachers Training College, Moscow, Part 2. p. 14–7. Lia O, Omre H, Tjelmeland H, Holden L, Egeland T. Uncertainties in reservoir production forecasts. Am Assoc Petrol Geol Bull 1997;81(5):775–802. Maslia ML, Prowell DC. Effect of faults on fluid flow and chloride contamination in a carbonate aquifer system. J Hydrol 1990;115: 1–49. Maas K. Drawdown pattern due to a well near a geological fault. In: Proceedings, 3rd international conference on the analytic element method in modeling groundwater flow, Brainerd, MN, USA, April 16–19, 2000.

1803

[35] Maas C. Groundwater flow to a well in a layered porous medium. 1. Steady flow. Water Resour Res 1987;23(8):1675–81. [36] Milne-Thomson LM. Hydrodynamical images. Proc Camb Philos Soc, Cambridge, England 1940;36:246–7. [37] Odling NE, Harris SD, Knipe RJ. Permeability scaling properties of fault damage zones in siliclastic rocks. J Struct Geol 2004;26:1727–47. [38] Pittman ED. Effect of fault-related granulation on porosity and permeability of quartz sandstones, Simpson Group (Ordovician), Oklahoma. AAPG Bull 1981;65(11):2381–7. [39] Shan C, Javendel I, Witherspoon PA. Characterization of leaky faults: study of water flow in aquifer–fault–aquifer systems. Water Resour Res 1995;31(12):2897–904. [40] Shipton ZK, Evans JP, Robeson KR, Forster CB, Snelgrove S. Structural heterogeneity and permeability in faulted eolian sandstone: implications for subsurface modeling of faults. AAPG Bull 2002;86(5):863–83. [41] Sigda JM, Goodwin LB, Mozley PS, Wilson JL. Permeability alteration in small-displacement faults in poorly lithified sediments; Rio Grande Rift, Central New Mexico. In: Haneberg WC et al., editors. Faults and subsurface fluid flow in the shallow crust. Geophysical monograph series, vol. 113. Washington, DC: AGU; 1999. p. 51–68. [42] Strack ODL. Groundwater mechanics. Englewood Cliffs, NJ: Prentice Hall; 1989. [43] Taylor WL, Pollard DD. Estimation of in situ permeability of deformation bands in porous sandstone, Valley of Fire, Nevada. Water Resour Res 2000;36(9):2595–606. [44] Yaxley LM. Effect of a partially communicating fault on transient pressure behavior. SPE Formation Eval 1987;4:590–8. [45] Zhang Y, Hobbs BE, Ord A, Barnicoat A, Zhao C, Walshe JL, et al. The influence of faulting on host rock permeability, fluid flow and ore genesis of gold deposits: a theoretical 2D numerical model. J Geochem Explorat 2003;78–79:279–84.