Two-dimensional seepage in porous media with heterogeneities

Journal of Geochemical Exploration 69–70 (2000) 251–255 www.elsevier.nl/locate/jgeoexp

Two-dimensional seepage in porous media with heterogeneities A.R. Kacimov a,*, Yu.V. Obnosov b a

Department of Soil and Water Sciences, Sultan Qaboos University, P.O. Box 34, Al-Khoud, Oman b Institute of Mathematics and Mechanics, Kazan University, Kazan, Russia

Abstract This paper deals with two-dimensional flow in heterogeneous media. New explicit analytical solutions are obtained for flow refraction on a growth fault, a growth anticline, and for a gravity–capillarity driven flow in an unsaturated medium with parabolic inclusion. Flow kinematic characteristics are studied. 䉷 2000 Elsevier Science B.V. All rights reserved. Keywords: conductivity; refraction; specific discharge; heterogeneity

1. Introduction Stratigraphic macroheterogeneities of permeability and porosity distributions are typical for geological formations and appear as growth faults, anticlines, diapirs, lenses, cavities, among others (Chapman, 1981, 1983). Microheterogeneities are less regular in shape and are by convention represented as random “patches” or pixels (“voxels”) with assigned step-wise constant physical characteristics (Barton and Beier, 1995; Torquato et al., 1999). Analytical solutions for flows through heterogeneous formations are rare because the refraction conditions along interfaces are difficult to satisfy in a rigorous way (Kacimov and Obnosov, 1997). Modeling of aquifer lenses dates back to the seminal ideas of Maxwell who showed that a uniform incident flow refracted by an ellipsoidal inclusion remains uniform within the inclusion with a three-dimensional (3D) distortion of streamlines in the ambient near-interface zone. Parquet-type heterogeneities composed of the so-called elementary cells repeating regularly in space are widely used in model* Corresponding author. E-mail addresses: [email protected] [email protected] (Y.V. Obnosov).

(A.R.

Kacimov),

ing subsurface flows (Renard and de Marsily, 1997). However, again, fine near-interface flow peculiarities are usually neglected. Obnosov (1996, 1999) developed a new approach to study refraction problems. The method enables one to determine 2D velocity fields in parquets, formations with single and multiple lenses at arbitrary conductivities of constituting porous components. Below we employ these solutions to analyze three cases. First, we study saturated flow near a “perfect” fault AB (Fig. 1) such that far from the fault, path lines are straight in intermittent layers while near AB they meander. Second, we consider saturated flow refracted by a parabolic inclusion (Fig. 4) with arbitrary orientation of the imposed gradient. Third, we develop the approach of Philip (1998) and study an unsaturated descendant flow near a parabolic lens (Fig. 6). 2. Saturated flows We consider a fault (Fig. 1a) composed of layers with constant thickness h and conductivities k1 and k2. The left-hand part of the massif of layers is shifted with respect to the right-hand part by a distance h along the fault line AB. We assume for definiteness

0375-6742/00/$ - see front matter 䉷 2000 Elsevier Science B.V. All rights reserved. PII: S0375-674 2(00)00131-X

252

A.R. Kacimov, Y.V. Obnosov / Journal of Geochemical Exploration 69–70 (2000) 251–255

Fig. 1. Refraction on a shift.

that k2 ⬎ k1 : Far from AB the specific discharge within each layer is constant. We specify its value V~ 1 …a; b† in the first medium. Then, in the second medium we have V~ 2 …k2 a=k1 ; b† at 兩x兩 ! ∞; i.e. a piece-wise 1D flow alternates its direction and magnitude when passing through the interfaces y ˆ 0; ^h; ^2h; … (Bear, 1972, p. 267; Miyazaki, 1993, pp. 81–87). At x ! 0 seepage becomes 2D, the shift distorts the straight trajectories. We originate a system of coordinates xOy at one of the corner points and orient Ox along the corresponding horizontal interface. We split out an elementary cell (Fig. 1b), which has a thickness of 2h and includes four layers extending to x ! ^∞: Within components h1 and h2, the hydraulic heads h1 and h2 are harmonic functions. Along the interfaces, the head and normal flux are continuous and hence, the conjugation

condition holds ~v1n …x; y† ˆ ~v2n …x; y†;

~v1t …x; y† v~ …x; y† ˆ 2t k1 k2

where t and n designate the tangential and normal components, respectively, of the specific discharge ~v1 and ~v2 in the two contacting media. The two Cartesian components of ~v1 …v1x ; v1y † are (Kacimov et al., 1999) v1x …x; y† ˆ

cos…2lV ⫹ g† …a 1 G ⫹ b1 =G† 1⫹D

sin…2lV ⫹ g† …a 1 G ⫺ b1 =G† v1y …x; y† ˆ ⫺ 1⫹D

…2†

where

Dˆ

k2 ⫺ k1 ; k1 ⫹ k2

p D 0 ˆ 1 ⫺ D2 ;

Aˆ

a1 ˆ …aD 0 ⫺ b†sin g; b1 ˆ …aD 0 ⫹ b†sin g; p sin g ˆ …1 ⫹ D†=2

Fig. 2. Streamlines in an elementary cell at k2 =k1 ˆ 5; a ˆ b ˆ 1:

…1†

Gˆ

p !2l sinh2 …px=h† ⫹ sin2 …py=h† cosh…px=h† ⫺ cos…py=h†

k1 ⫹ k2 2k2

A.R. Kacimov, Y.V. Obnosov / Journal of Geochemical Exploration 69–70 (2000) 251–255

Fig. 3. Nondimensional travel time through an elementary cell as a function of the initial location of a marked particle.

l ˆ p⫺1 arcsin D sin…py=h† V ˆ ⫿arcsin p 2 sinh …px=h† ⫹ sin2 …py=h† and “ ⫺ ” corresponds to x ⬎ 0 and “ ⫹ ” corresponds to x ⬍ 0: The components of ~v2 are similar to Eq. (2). Using Mathematica (Wolfram, 1991) we integrated Eq. (1) as a system of two ordinary differential equations. Fig. 2 shows the streamlines within one elementary cell 0 ⬍ y ⬍ 2h at a ˆ b ˆ 1 and k2 =k1 ˆ 5: Lines 1 and 5 located far from the shift (clearly, at the macroscale in Fig. 1a all streamlines eventually intersect the shift) behave “normally”, i.e. are not

253

curved. Streamline 3 is nearly not influenced by the layering, it passes from one strip of high conductivity into another without meandering as it would be in a homogeneous medium of conductivity k2. Streamlines 2 and 4 have nontrivial shapes because the corresponding particles are most affected by the vertical and horizontal interfaces. Fig. 3 illustrates a nondimensional travel time k1 t=h as a function of x=h where t is the time of passage of a marked particle through an elementary cell (additionally, we selected the porosity of both layers to be m ˆ 0:3†: The extrema of this curve are similar to Kacimov et al. (1999). In particular, the minimum corresponds to particle 3 in Fig. 2, which is “funneled” by the shift. Consider now a parabolic lens with an interface y2 ˆ l…2x ⫹ l† separating two zones of constant conductivity k1 and k2 outside and inside the lens, respectively (Fig. 4). The specific discharge in the ~ x ; Vy †; first medium far away from the interface, V…V is given. The orientation of the incident flow, a i can be arbitrary. In particular, the parabola limbs can be turned upward that corresponds to a common layering in sedimentary basins. Along the interface we have the same refraction condition (1). The goal is to determine two vector fields ~v…v1x ; v1y † and ~v2 …v2x ; v2y †; which depend on k ˆ k2 =k1 ; a i and l. We introduce a complex coordinate z ˆ x ⫹ iy and reformulate condition (1) using the parabola arc coordinate. After some algebra, the components of the

Fig. 4. Parabolic lens in saturated flow.

254

A.R. Kacimov, Y.V. Obnosov / Journal of Geochemical Exploration 69–70 (2000) 251–255

How quickly will the particle reach the parabola tip? Since v1y ˆ 0 the time of advection for this particle is Z⫺ l=2 dx …4† tˆm ⫺ d v1x …x†

Fig. 5. Nondimensional travel time along the parabola axis at Vy ˆ 0 as a function of the abscissa of the starting point.

specific discharge v1x …x; y†; v1y …x; y†; v2x …x; y† and v2y …x; y† are p k ⫺ 1 p p …k r ⫺ xVx ⫺ sign‰yŠ r ⫹ xVy † v1x ˆ Vx ⫹ l 2rk p k ⫺ 1 p p …sign‰yŠk r ⫹ xVx ⫹ r ⫺ xVy † v1y ˆ Vy ⫺ l 2rk v2x ˆ kVx ;

v2y ˆ Vy

…3† p where r ˆ 兩z兩 ˆ x2 ⫹ y2 : As Eq. (3) shows, the velocity in the lens is constant but its magnitude and direction differ from the incident velocity as was the case for the Maxwell ellipsoids. Eq. (3) can be used to analyze the kinematic characteristics similarly to the shift. For example, let us consider Vy ˆ 0 and a marked particle moving along the parabola axis from point D…x0 ˆ ⫺d† towards the lens (Fig. 4).

Fig. 6. Descending unsaturated flow near a parabolic boundary.

we can easily integrate Eq. (4) (the corresponding explicit expression is dropped for brevity). Fig. 5 shows a nondimensional travel time k1 t=…lm† as a function of d=l at k ˆ 0:1; 0.5 and 5 (curves 1–3, respectively). Analogously, from Eq. (3) we can assess the value of dispersivity, which in the advective dispersion model is proportional to 兩~v兩: Eq. (4) shows that the lens accelerates at k ⬎ 1 and decelerates at k ⬍ 1 marked particles proportionally to k ⫺ 1 and, hence, the nonuniformity of the Darcian velocity in the ambient medium is significant even for moderate conductivity contrasts. On the other hand, dispersivity within the parabola is constant and migration patterns there can be reconstructed from the known 2D transport equation with constant coefficients (Bear, 1972). For unsaturated conditions, Philip (1998) solved the refraction problem for a parabolic lens, which axis is cooriented with a descending flow. He showed that that flow is uniform in the interior zone. Hence, the parabola is a unique interface because for ellipsoids and ellipses an unsaturated flow cannot be uniform inside the lens. 3. Unsaturated flow We consider the unsaturated flow to a parabolic boundary G . Unlike the previous problem, the axis is strictly vertical and the equation of the parabola is 2lZ ˆ X 2 ⫺ l2 : Sufficiently far from G , the flow is uniform with a constant velocity U∞ oriented downward (Fig. 6). Philip et al. (1989) discussed this pattern in application to subterranean holes, tunnels or caves disturbing natural infiltration. Philip (1998) illustrated how an isobaric or impermeable parabolic cavity can be interpreted as a capillary–barrier interface. We apply a new technique to the Philip et al. (1989) model and derive a general solution, which is applicable to numerous flow conditions. According to the so-called quasilinear model (Philip, 1998) the relation between the moisture potential C and unsaturated conductivity k is k…C† ˆ k s exp…aC†

A.R. Kacimov, Y.V. Obnosov / Journal of Geochemical Exploration 69–70 (2000) 251–255

where ks is the saturated conductivity and a is the sorptivity length. The Kirchhoff potential Q…X; Z† ˆ k…C†=a satisfies the advective dispersion equation. Along G we specify an arbitrary distribution of Q . We introduce the parabolic coordinates …z; h† through X ˆ zh; Z ˆ …z 2 ⫺ h2 †=2 in which the governing equation becomes   22 Q 22 Q 2Q 2Q ⫹ ⫺ 2s z ⫺h ˆ0 …5† 2z 2h 2h 2 2z 2 where s ˆ 0:5al and uniform infiltration is specified by a given value of Q ∞ ˆ U∞ =a at z ! ∞: Philip et al. (1989) found a particular solution to Eq. (5) for Q ˆ Q i ˆ const along G . If Q G varies, then a general solution can be derived in terms of the Hermite polynomials Hn and an associated function In. For example, if

Q G ˆ Q i ⫺ b exp…⫺g2 z2 † where b and g are given constants, then the potential distribution is

Q ˆ Q ∞ ⫹ …Q i ⫺ Q ∞ † ⫺b

∞ X

En …h† ˆ

…⫺1†n g2n H2n …z†I2n …h† p I2n … s† 2 n!…1 ⫹ g2 †n⫹0:5 2n

nˆ0

where

erfc h p erfc s

(

Hn …⫺ih†

…6†

for even n

⫺iHn …⫺ih† for odd n

I n … h† ˆ E n … h†

Z∞ h

e⫺t En⫺2 …t† dt 2

At b ˆ 0 Eq. (6) reduces to the Philip (1998) solution. 4. Conclusions The solutions derived give a better insight of the 2D fluid movement in structurally nonuniform porous media. Eqs. (2) and (3) for the specific discharge in saturated flow and Eq. (6) for the Kirchhoff potential in unsaturated quasilinear flow can be used in solute

255

transport models and other algorithms of particle tracking (e.g. the random walk method). These solutions circumvent mesh generation and numerical differentiation (with associated errors) of head (pressure) values as in standard FDM–FEM codes. Acknowledgements This study was supported by Sultan Qaboos University, project AGSWAT 9903 and by Russian Foundation of Basic Research, grants N99-0100364, 99-01-00173. Helpful comments of an anonymous referee are appreciated. References Barton, H.H., Beier, R.A., 1995. Fractals in Petroleum Geology and Earth Processes, Plenum Press, New York. Bear, J., 1972. Dynamics of Fluids in Porous Media, Elsevier, New York. Chapman, R.E., 1981. Geology and Water, Nijhoff, The Hague. Chapman, R.E., 1983. Petroleum Geology, Elsevier, Amsterdam. Kacimov, A.R., Obnosov, Yu.V., 1997. Explicit solution to a problem of seepage in a checker-board massif. Transport Porous Media 28, 109–124. Kacimov, A.R., Obnosov, Yu.V., Yakimov, N.D., 1999. Groundwater flow in a medium with a parquet-type conductivity distribution. J. Hydrol. 226, 242–249. Miyazaki, T., 1993. Water Flow in Soils, Marcel Dekker, New York. Obnosov, Yu.V., 1996. Exact solution of a problem of R-linear conjugation for a rectangular checkerboard field. Proc. R. Soc. London A 452, 2423–2442. Obnosov, Yu.V., 1999. Periodic heterogeneous structures: new explicit solutions and effective characteristics of refraction of an imposed field. SIAM J. Appl. Math. 59, 1267–1287. Philip, J.R., 1998. Seepage shedding by parabolic capillary barriers and cavities. Water Resour. Res. 34, 2827–2835. Philip, J.R., Knight, J.H., Waechter, R.T., 1989. The seepage exclusion problem for a parabolic and paraboloidal cavities. Water Resour. Res. 25, 605–618. Renard, Ph., de Marsily, G., 1997. Calculating equivalent permeability: a review. Adv. Water Res. 20, 253–278. Torquato, S., Kim, I.C., Cule, D., 1999. Effective conductivity, dielectric constant, and diffusion coefficient of digitized composite media via first-passage-time equations. J. Appl. Phys. 85, 1560–1571. Wolfram, S., 1991. Mathematica. A System for doing Mathematics by Computer, Addison-Wesley, Redwood City, CA.