Convergence to the macroscopic model is shown by extracting the weak limits of the microscopic model solutions

Electronic Journal of Differential Equations, Vol. 2010(2010), No. 94, pp. 1–19.

ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu

A MODEL FOR SINGLE-PHASE FLOW IN LAYERED POROUS MEDIA

DANIEL J. COFFIELD JR., ANNA MARIA SPAGNUOLO

Abstract. Homogenization techniques are used to derive a double porosity model for single phase flow in a reservoir with a preferred direction of fracture.

The equations in the microscopic model are the usual ones derived from Darcy’s law in the fractures and matrix (rock). The permeability coefficients over the matrix domain are scaled, using a parameter ε, based on the fracture direction in the reservoir. The parameterε represents the size of the parts of the matrix blocks that are being homogenized and the scaling preserves the physics of the flow between matrix and fracture as the blocks shrink.

Convergence to the macroscopic model is shown by extracting the weak limits of the microscopic model solutions. The limit (macroscopic) model consists of Darcy flow equations in the matrix blocks and fracture sheet, with additional terms in the fracture sheet equation. Together, these terms represent the fluid exchange between the matrix blocks and the fracture sheet.

1. Introduction

It is well known that fluid flows in fractured reservoirs as if the reservoir has two porous structures, one associated with the fissure system and the other as- sociated with the porous rocks. From this came the dual (double)-porosity (dual- permeability) concept. Models of this type for single-phase flow were first developed in [4, 13, 14, 19] using physical arguments. Since then, a more general form of a double-porosity model for single-phase flow in totally fractured reservoirs was de- rived on physical grounds and by using formal homogenization in [2, 3, 7, 10], and then proven rigorously using homogenization techniques [1]. For the main ideas of fluid flow in porous media and homogenization theory, see [5, 15], respectively. In addition, a general treatment of homogenization and porous media can be found in [12].

Double-porosity models have also been used to better model partially fissured media [9]. Additionally, the dual-porosity concept has been applied to layered porous media. In [8] a dual-porosity model for two-phase immiscible flow in a vertically fractured reservoir is derived using formal homogenization techniques.

A model of dual-porosity type is derived on physical grounds in [6] for a porous

2000Mathematics Subject Classification. 76S05, 35B27.

Key words and phrases. Layered media; naturally-fractured media;

double-porosity model; homogenization.

c

2010 Texas State University - San Marcos.

Submitted February 23, 2010. Published July 13, 2010.

1

medium with relatively thin rock cells. Furthermore, multiple-porosity models have been derived using homogenization techniques in [1, 11, 16, 17, 18].

Fractured reservoirs are often classified according to the extent and character- istics of the fracture system. In a totally fractured porous medium, the fracture system is fully developed and forms a connected set. The fractures divide the rocks (matrix blocks) into disconnected pieces so that for fluid to flow from one matrix block to another, it must first pass through the fracture system. In some sense, a totally fractured reservoir is an idealization because most reservoirs will have a less-developed fracture system that is not completely connected. Reservoirs of this type are called partially fractured (or partially fissured) porous media. The matrix blocks are not necessarily disconnected in this case and the fractures are not nec- essarily connected. A special case of fractured reservoirs occurs when the fractures are naturally developed with a preferred direction. Often the fractures occur in planes with a particular orientation so that the matrix blocks, separated by the fractures, form a layered medium.

In this work, we consider single-phase flow through a totally fractured layered medium. The porous medium consists of horizontal fractures so that the elongated matrix blocks are stacked vertically. This is a typical structure found in certain sedimentary rocks such as shale. The primary contribution of this work is the derivation of a well-posed dual-porosity model for flow through such a reservoir by using standard fluid flow equations posed at the micro-scale and then averaging this flow through rigorous homogenization techniques. We followed the techniques in [1], but, in our case, the homogenization is not done uniformly. More specifically, we homogenize in the vertical direction only, reflecting the physical nature of the reservoir: the structure is periodic in only one (the vertical) direction. The one- dimensional homogenization technique also requires careful component-wise (not uniform) scaling of the matrix permeability. Scaling of this type has been used in [8]. Unique to our problem is that vertical side fractures do not contain any rock and therefore are not averaged, but they are coupled to the flow equations that are being homogenized. Also, a unique feature of our final homogenized system is that horizontal flow occurs simultaneously in the fracture sheet and matrix blocks.

In what follows, we describe the fractured porous medium and the assumptions in Section 2. In Section 3, we describe the microscopic model. Section 4 includes many technical lemmas and Section 5 defines the macroscopic limit model and includes the details of how the microscopic model converges to the limit model.

Finally, Section 6 briefly discusses some extensions to this work.

2. Preliminaries

We consider the reservoir Ω = Ω¹×Ω²×Ω³⊂R³to be a bounded, rectangular domain that is the union of the closures of congruent rectangular domains. That is, we denote the standard rectangular cell byQand Ω is the interior of the union of translations ofQ:

Ω =∪_α∈I(1)(Q+cα),

where thecα’s are translations and I(1) ={0,1,2, . . . , N(1)} is an index set (See Fig. 1).

Additionally, we require the cells to be disjoint: (Q+cα)∩(Q+cβ) = ∅ for α6=β. LetC(1) ={cα:α∈I(1)}be the set of translations. Then (0,0,0)∈C(1) and we call itc₀. Also, there is a c_α= (0,0, r)∈C(1), call itc₁, withr∈R, such

Figure 1. The reservoir Ω = Ω¹×Ω²×Ω³.

thatC(1) ={αc1:α∈I(1)}. That is,C(1) is made up of integer multiples of the standard translation,c1, and the translations only occur in the third component.

The standard rectangular cell Q =Q¹×Q²×Q³ consists of three parts: the matrix block, denoted byQ_m; the fracture surrounding the matrix, denoted byQ_f; and the internal boundary betweenQ_mandQ_f, denoted by∂Q_m(See Fig. 2). The domainQ_mis itself rectangular,Q_m=Q¹_m×Q²_m×Q³_m, and is compactly contained in Q. The domain Q_f is connected and ∂Q_m is Lipschitz and piecewise smooth.

We also define Q¹_f =Q¹\Q¹_m, Q²_f =Q²\Q²_m, and Q³_f =Q³\Q³_m. The domainQ³_f consists of two pieces, which we denote byQ^3,1_f and Q^3,2_f , with |Q^3,1_f |=|Q^3,2_f |=

1 2|Q³_f|.

Figure 2. The standard cellQand its subsets.

To homogenize, we will shrink the cells in the direction of periodicity, the vertical direction. We introduce the notation,^εQ=Q¹×Q²×εQ³, and call this theε-cell.

Similarly, we define^εQm=Q¹_m×Q²_m×εQ³_mand^εQf ={(y1, y2, εy3) : (y1, y2, y3)∈ Qf}. When ε= 1, we have ^εQ= Q, the standard cell. As ε tends to 0 the cells shrink linearly in the vertical direction only. For each ε, we have a new index set I(ε) ={0,1,2, . . . , N(ε)}and set of translationsC(ε) ={cα:α∈I(ε)}so that the ε-cells fill up Ω (See Fig. 3).

Figure 3. Vertical cross-sections through the middle of Ω showing the vertical homogenization asε→0 and theε-dependent subsets of Ω.

Next, we define a matrix domain and fracture domain for eachε >0 as follows:

Ω^ε_m= Ω∩ ∪_α∈I(ε)(^εQ_m+εc_α) and Ω^ε_f = Ω∩ ∪_α∈I(ε)(^εQ_f+εc_α), (See Fig. 3). Note that the translations are 0 in the first two components so theε only affects the third component of the translationcα. For convenience we assume that the sequence ofε’s are such that∂Ω⊂∂Ω^ε_f.

Next, we introduce notation for various components of Ω. Since, Ω¹=Q¹ and Ω² = Q², we similarly define the following: Ω¹_m = Q¹_m, Ω¹_f = Q¹_f, Ω²_m = Q²_m, and Ω²_f =Q²_f. For simplicity, we define Ω^1,2 = Ω¹×Ω², Ω^1,2_m = Ω¹_m×Ω²_m, and Ω^1,2_f = (Ω¹_f×Ω²)∪(Ω¹×Ω²_f). To describe the setting of the macroscopic model, we decompose Ω into three pieces: the fracture sheet, denoted by ΩM = Ω^1,2_m ×Ω³, corresponding to the part of Ω with matrix blocks in it; the side fractures, denoted by ΩF = Ω^1,2_f ×Ω³, that do not shrink asεtends to 0 and correspond to the fractures on the lateral sides of the fracture sheet; and the internal boundary between the two, denoted by∂Ω^1,2_m ×Ω³. Note that Ω^1,2_m ×Ω^3,ε_m = Ω^ε_mand Ω^1,2_m ×Ω^3,ε_f = Ω^ε_f\ΩF, where Ω^3,ε_m and Ω^3,ε_f are twoε-dependent subsets of Ω³(See Fig. 4).

As the cells shrink, their horizontal components remain unchanged. That is,Q¹, Q², and their subsets are not shrinking at all. So, while the vertical components of Ω andQare on different scales asεtends to 0, the horizontal components remain on the same scale and are, in fact, the same: Qⁱ = Ωⁱ, Qⁱ_m= Ωⁱ_m, andQⁱ_f = Ωⁱ_f, for i = 1,2. Because of this, we ignore the redundant spaces for the horizontal components and refer only to Ω¹, Ω², and their subsets.

For each 0< ε≤1, we define a functionc^ε: Ω³→C(ε) by letting c^ε(x3) be the translation corresponding to theε-cell that contains Ω^1,2× {x3}. The function c^ε can be properly defined by assigning the ε-cell boundaries to either the cell above or below it in a non-overlapping way. Finally, we let [0, T] be the time interval of interest and useJ to denote the open time interval (0, T).

Figure 4. A vertical cross-section through the middle of Ω show- ing theε-dependent subsets of Ω³.

3. The Microscopic Model

In this section we pose the model on the microscopic scale. Forx= (x1, x2, x3), let ρ^ε(x, t) and σ^ε(x, t) be the fluid density in Ω^ε_f and Ω^ε_m respectively. Next, we assume the fluid is a liquid with Newtonian viscosity, µ, and that it satisfies the following equations of state: dρ^ε =cρ^εdp^ε and dσ^ε = cσ^εdp^ε, where c is the constant (small) compressibility andp^εis the pressure.

Letφ(x) andk(x) be the porosity and permeability matrix, respectively, in the blocks of the original reservoir, Ω^ε=1_m . Since the blocks are periodically distributed over Ω in the vertical direction, we assume thatφandkare vertically periodic with a period ofQ³. We then letφ^ε(x) =φ(x₁, x₂,^x_ε³) andk^ε(x) =k(x₁, x₂,^x_ε³) be the porosity and permeability matrix in Ω^ε_mso thatφ^εandk^εareεQ³-periodic in the vertical direction.

Following [1], the fluid flow is assumed to be governed by Darcy’s law in Ω^ε_mand Ω^ε_f. We define Φ^∗(x) and K^∗(x) to be the porosity and scalar permeability of the original fracture domain, defining them first over all of Ω and then considering their restrictions to Ω^ε=1_f . Now, φ, Φ^∗, k andK^∗ are all smooth and bounded, with φ, Φ^∗, andK^∗being uniformly positive. Also,kis uniformly bounded component-wise and uniformly positive definite.

Now we describe the microscopic model. Since we are only homogenizing in the one (vertical) direction of periodicity, we must carefully scale the permeability matrixk^ε. We use a scaling similar to the one done in [8]. LetH be the matrix

H=





1 0 0 0 1 0 0 0 ε



.

We then define a scaled permeability matrix,bk^ε, bybk^ε=Hk^εH. So,

bk^ε(x) =





k^ε₁₁(x) k^ε₁₂(x) εk₁₃^ε (x) k^ε₂₁(x) k^ε₂₂(x) εk₂₃^ε (x) εk^ε₃₁(x) εk^ε₃₂(x) ε²k^ε₃₃(x)





=





k11(x1, x2,^x_ε³) k12(x1, x2,^x_ε³) εk13(x1, x2,^x_ε³) k₂₁(x₁, x₂,^x_ε³) k₂₂(x₁, x₂,^x_ε³) εk₂₃(x₁, x₂,^x_ε³) εk₃₁(x₁, x₂,^x_ε³) εk₃₂(x₁, x₂,^x_ε³) ε²k₃₃(x₁, x₂,^x_ε³)



.

Throughout we use ν to mean the outward pointing unit normal. Then using the conservation of mass equations combined with Darcy’s law and the equations of state, the fluid flow in Ω^ε_f and Ω^ε_mcan be described in the standard way as follows:

Φ^∗ρ^ε_t− ∇ · K^∗

µc∇ρ^ε

=f in Ω^ε_f, t >0, (3.1) K^∗

µc∇ρ^ε·ν = bk^ε

µc∇σ^ε·ν on∂Ω^ε_m, t >0, (3.2) K^∗

µc∇ρ^ε·ν = 0 on∂Ω, t >0, (3.3) ρ^ε=ρ_init in Ω^ε_f, t= 0, (3.4) φ^εσ^ε_t− ∇ · bk^ε

µc∇σ^ε

!

=f in Ω^ε_m, t >0, (3.5) σ^ε=ρ^ε on∂Ω^ε_m, t >0, (3.6) σ^ε=ρinit in Ω^ε_m, t= 0. (3.7) The function f = f(x, t) represents external sources and sinks. The boundary conditions (3.2) and (3.6) respectively represent conservation of mass and continuity of pressure between Ω^ε_f and Ω^ε_m. Also, (3.3) gives no-flow across the external boundary ∂Ω. For the sake of convenience, we are ignoring gravity in the model but could easily account for it.

The need for the ε-scaling in bk^ε can be seen by considering matrix-to-fracture flow as the blocks shrink. The scaling causes the matrix blocks to be less permeable asε→0, and if no scaling is done, the matrix-to-fracture flow will blow up as the blocks shrink. Using a similar technique to that found in [2], we choose theε-scaling in bk^ε to conserve matrix-to-fracture flow in a certain sense. The original (ε = 1) matrix-to-fracture flow is given by:

Z

∂Ω^ε=1_m

k(x)

µc ∇σ^ε=1(x)·ν ds(x) = Z

Ω^ε=1_m

∇ · k(x)

µc ∇σ^ε=1(x)

dx

=

N(1)

X

i=0

Z

Q_m+c_i

∇ · k(x)

µc ∇σ^ε=1(x)

dx,

(3.8)

whereN(1) is the number of cells in the original (ε= 1) domain. Now, for arbitrary ε <1, we consider the matrix-to-fracture flow without scaling. Making a change of

variables and recalling thatkisQ³-periodic, the flow is given by Z

∂Ω^ε_m

k^ε(x)

µc ∇σ^ε(x)·ν ds(x)

= Z

Ω^ε_m

∇ ·

k^ε(x) µc ∇σ^ε(x)

dx

=

N(ε)

X

i=0

Z

εQ_m+εc_i

∇ ·

k^ε(x) µc ∇σ^ε(x)

dx

=ε

N(ε)

X

i=0

Z

Qm+ci

∇^1/ε_z ·

k(z1, z2, z3)

µc ∇^1/ε_z σ^ε(z1, z2, εz3)

dz,

(3.9)

where ∇^1/εz means _∂z^∂

1,_∂z^∂

2,¹_ε∂z^∂

3

. Comparing (3.8) and (3.9), we see that the extraεis accounted for sinceN(ε) is on the order ofε⁻¹N(1). But, the extra ¹_ε’s are not accounted for. By scalingk^ε as we have done inbk^ε, we account for them and keep the flow from blowing up. Thus, choosing bk^ε =Hk^εH seems to be an appropriate choice for conserving the matrix-to-fracture flow as we homogenize.

To simplify things, we introduce more notation:

Φ(x) =

(Φ^∗(x) ifx∈ΩF

|Q³_f|

|Q³|Φ^∗(x) ifx∈ΩM, Λ(x) = ΛF(x) =K^∗(x) µc ,

Λ_M(x) = |Q³_f|

|Q³|

K_M(x)

µc , λ^ε(x) = k^ε(x)

µc , bλ^ε(x) = bk^ε(x) µc .

Throughout we assume that f ∈ L²(J;L²(Ω)), ρ_init ∈ L²(Ω), and use (·,·)_X to denote the standard inner product of L²(X). Then for ϕ∈L²(J;H¹(Ω)), we use the boundary conditions and (3.5) to see that

−(∇ ·(Λ∇ρ^ε), ϕ)Ω^ε_f = (Λ∇ρ^ε,∇ϕ)Ω^ε_f −(Λ∇ρ^ε·νf, ϕ)∂Ω^ε_m−(Λ∇ρ^ε·ν, ϕ)∂Ω

= (Λ∇ρ^ε,∇ϕ)Ω^ε_f + (bλ^ε∇σ^ε·νm, ϕ)∂Ω^ε_m

= (Λ∇ρ^ε,∇ϕ)Ω^ε_f + (φ^εσ^ε_t, ϕ)Ω^ε_m−(f, ϕ)Ω^ε_m+ (bλ^ε∇σ^ε,∇ϕ)Ω^ε_m, where ν_f and ν_m are the outward pointing unit normals from Ω^ε_f and Ω^ε_m respec- tively. So a weak form for the microscopic model is

(Φ^∗ρ^ε_t, ϕ)Ω^ε_f×J+ (Λ∇ρ^ε,∇ϕ)Ω^ε_f×J+ (φ^εσ_t^ε, ϕ)Ω^ε_m×J+ (bλ^ε∇σ^ε,∇ϕ)Ω^ε_m×J

= (f, ϕ)_Ω×J, ϕ∈L²(J;H¹(Ω)),

(3.10) (φ^εσ^ε_t, ψ)Ω^ε_m×J+ (bλ^ε∇σ^ε,∇ψ)Ω^ε_m×J= (f, ψ)Ω^ε_m×J, ψ∈L²(J;H₀¹(Ω^ε_m)), (3.11) ρ^ε=σ^ε forx∈∂Ω^ε_m, t >0, (3.12) ρ^ε=ρ_init, σ^ε=ρ_init forx∈Ω, t= 0. (3.13) Theorem 3.1. For each ε, there exists a unique solution, ρ^ε and σ^ε, to the sys- tem (3.10)–(3.13). Additionally, ρ^ε ∈ H¹(J;L²(Ω^ε_f))∩L²(J;H¹(Ω^ε_f)) and σ^ε ∈ H¹(J;L²(Ω^ε_m))∩L²(J;H¹(Ω^ε_m)).

Proof. First, letχ^ε_f be the characteristic function of Ω^ε_f and θ^ε(x) =

(ρ^ε forx∈Ω^ε_f σ^ε forx∈Ω^ε_m,

α^ε=χ^ε_fΦ^∗+ (1−χ^ε_f)φ^ε, κ^ε=χ^ε_fΛ + (1−χ^ε_f)bλ^ε. Then (3.10) and (3.13) can be rewritten as

(α^εθ_t^ε, ϕ)Ω×J+ (κ^ε∇θ^ε,∇ϕ)Ω×J = (f, ϕ)Ω×J, ∀ ϕ∈L²(J;H¹(Ω)), (3.14) θ^ε=ρinit forx∈Ω, t= 0. (3.15) This is a single parabolic problem with discontinuous coefficients. It is well known that this problem has a unique solution in H¹(J;L²(Ω))∩L²(J;H¹(Ω)). By re- stricting the solution to Ω^ε_f and Ω^ε_m, we obtainρ^εandσ^ε respectively.

Hereafter we use C to denote a positive constant that is independent of ε, though not necessarily the same one at each occurrence. Also, we use∇^εto denote (_∂x^∂

1,_∂x^∂

2, ε_∂x^∂

3). Then, from the standard energy estimates used to prove Theorem 3.1, we obtain the following result.

Lemma 3.2.

kρ^εkL²(J;L²(Ω^ε_f))+kσ^εkL²(J;L²(Ω^ε_m))≤C, kρ^ε_tk_L2(J;L²(Ω^ε_f))+kσ_t^εk_L2(J;L²(Ω^ε_m))≤C, k∇ρ^εk_L2(J;L²(Ω^ε_f))+k∇^εσ^εk_L2(J;L²(Ω^ε_m)) ≤C.

4. Technical Lemmas

To prove the main result we will need to use various technical lemmas which we state in this section. First, for a measurable function φon Ω, a subset of Ω, or a subset of the boundary of Ω, we define

φ(xe 1, x2, x3, y3) =φ(x1, x2, εy3+c^ε(x3)).

We call the operator mappingφtoφethe dilation operator. Recalling the definition of c^ε, we are actually defining a family of operators, one for each ε, 0 < ε ≤ 1.

However, for convenience we will only write “

e”, and leave the ε-dependence as implicitly understood. We first note that the dilation operator commutes with ad- dition and multiplication: φψf =φeψeandφ^+ψ=φ+e ψ. In the following lemmas, wee prove some deeper properties of the dilation operator that will be used extensively.

Lemma 4.1. If ψ∈H¹(Ω), then _∂y^∂

3ψe=ε]^∂

∂x₃ψ.

Lemma 4.2. If ψ, ϕ∈L²(Ω) (orH¹(Ω) when appropriate), then for r=m, f,or blank,

(ψ,e ϕ)e_Ω1,2

p(r)×Ω³×Q³_r =|Q³|(ψ, ϕ)_Ω^1,2

p(r)×Ω^3,εr , kψke _L2(Ω^1,2_p(r)×Ω³×Q³_r)=|Q³|^1/2kψk_L2(Ω^1,2_p(r)×Ω^3,εr ), k ∂

∂y3

ψke _L2(Ω^1,2_p(r)×Ω³×Q³_r)=ε|Q³|^1/2k ∂

∂x3

ψk_L2(Ω^1,2_p(r)×Ω^3,εr ),

(ψ, ϕ)e _Ω1,2

r ×Ω³×Q³ = (ψ,ϕ)e_Ω1,2

r ×Ω³×Q³,

whereΩ^3,ε means Ω³ andp(r) =

(m ifr=f or m, f ifr=blank.

Proof. Lettingz3=εy3+c^ε(x3), the first result is a calculation.

(ψ,e ϕ)e_Ω1,2

p(r)×Ω³×Q³_r = Z

Ω^1,2_p(r)

Z

Ω³

Z

Q³_r

ψ(x₁, x₂, εy₃+c^ε(x₃))ϕ(x₁, εy₃+c^ε(x₃))dy₃dx

= Z

Ω^1,2_p(r)

Z

Ω³

Z

εQ³_r+c^ε(x3)

ψ(x₁, x₂, z₃)ϕ(x₁, x₂, z₃)ε⁻¹dz₃dx

=|Q³| Z

Ω^1,2_p(r)

Z

Ω^3,εr

ψ(x1, x2, z3)ϕ(x1, x2, z3)dz3dx1dx2

=|Q³|(ψ, ϕ)_Ω^1,2

p(r)×Ω^3,εr .

The next two results follow from the first result and the previous lemma. The last result is also a calculation.

(ψ, ϕ)e _Ω1,2 r ×Ω³×Q³

= Z

Ω^1,2_r

Z

Ω³

Z

Q³

ψ(x1, x2, εy3+c^ε(x3))ϕ(x1, x2, x3)dy3dx

= Z

Ω^1,2r

Z

Ω³

Z

εQ³+c^ε(x3)

ψ(x1, x2, z3)ϕ(x1, x2, x3)ε⁻¹dz3dx

=ε⁻¹ Z

Ω^1,2_r

Z

Ω³

ψ(x1, x2, z3)Z

εQ³+c^ε(z₃)

ϕ(x1, x2, x3)dx3

dz3dx1dx2

=ε⁻¹ Z

Ω^1,2_r

Z

Ω³

ψ(x1, x2, z3)Z

Q³

ϕ(x1, x2, εy3+c^ε(z3))ε dy3

dz3dx1dx2

= (ψ,ϕ)e _Ω1,2

r ×Ω³×Q³.

Lemma 4.3. Lets=F orM andr=m,f, orblank. Ifϕ∈L²(Ωs)is considered as a function inL²(Ω_s×Q³_r)that is constant iny₃, thenϕe→ϕstrongly inL²(Ω_s× Q³_r)asε→0.

Proof. Using the Dominated Convergence Theorem, the result is straightforward forψ∈C₀^∞(Ωs):

ε→0lim Z

Ω_s×Q²_r

(ψ(x, ye ₂)−ψ(x))²dx dy₂

= Z

Ω_s×Q²_r

( lim

ε→0ψ(x₁, εy₂+c^ε(x₂))−ψ(x))²dx dy₂= 0.

The full result then follows from Lemma 4.2 and the fact that functions inC₀^∞(Ωs)

are dense inL²(Ωs).

Corollary 4.4. We have the following estimates:

kfσ^εk_L2(Ω³;H¹(Ω^1,2m×Q³_m×J))≤C, kρe^εk_L2(Ω³;H¹(Ω^1,2_f ×Q³×J))≤C, kρe^εk_L2(Ω³;H¹(Ω^1,2_m×Q³_f×J)) ≤C,

k ∂

∂y₃ρe^εk_L2(Ω_F;L²(Q³;L²(J)))≤εC, k ∂

∂y3

ρe^εk_L2(Ω_M;L²(Q³_f;L²(J)))≤εC.

The proofs for the above estimates are straighforward calculations using Lemmas 3.2 and 4.2.

5. The Main Result

When we homogenize, the microscopic model converges to an averaged, limit model. We now describe that limit model which we call the macroscopic model.

Figure 5. The macroscopic model domains, ΩF and ΩM, includ- ing a horizontal cross-section.

Let ρ^F(x, t) and ρ^M(x, t) be the macroscopic, averaged fluid density in side fractures and fracture sheet, Ω_F and Ω_M, respectively (See Fig. 5). Also, let

ρ(x, t) =

(ρ^F(x, t) ifx∈ΩF,

ρ^M(x, t) ifx∈ΩM, Λ^#(x) =

(ΛF(x) ifx∈ΩF, ΛM(x) ifx∈ΩM. We let∇= (_∂x^∂

1,_∂x^∂

2,_∂x^∂

3),∇x₁,x₂= (_∂x^∂

1,_∂x^∂

2), and∇x₁,x₂,y₃ = (_∂x^∂

1,_∂x^∂

2,_∂y^∂

3).

Then, withf_m,f_M¹ , andf_M² defined as below, ρ^F andρ^M satisfy the following:

Φρ^F_t − ∇ · ΛF∇ρ^F

=f in ΩF, t >0, (5.1) Φρ^M_t − ∇x1,x2· Λ_M∇x1,x2ρ^M

− ∂

∂x1

f_M¹ − ∂

∂x2

f_M² =f+f_m in Ω_M, t >0, (5.2) ΛF∇ρ^F·ν= 0 on∂Ω\(Ω^1,2_m ×∂Ω³), t >0, (5.3)

ρ=ρ_init in Ω, t= 0, (5.4)

∂

∂x₁ρ^F =|Q³_f|

|Q³|

∂

∂x₁ρ^M, f_M¹ = 0 on∂Ω¹_m×Ω²_m×Ω³, t >0, (5.5)

∂

∂x₂ρ^F =|Q³_f|

|Q³|

∂

∂x₂ρ^M, f_M² = 0 on Ω¹_m×∂Ω²_m×Ω³, t >0. (5.6) The boundary condition (5.3) holds on the boundary of Ω, excluding the top and bottom of ΩM. The conditions (5.5) and (5.6) hold on the internal boundary between ΩF and ΩM. Asε tends to zero, the blocks shrink to become horizontal plates. There are infinitely many of them, one for each x3. Let σ(x, y3, t) be the macroscopic fluid density in the matrix blocks and note that the macroscopic poros- ity and scaled permeability matrix are φ(x₁, x₂, y₃) and λ(x₁, x₂, y₃) respectively.

Thenσsatisfies the following:

φσt− ∇x₁,x₂,y₃·(λ∇x₁,x₂,y₃σ) =f(x, t) in ΩM×Q³_m, t >0, (5.7) σ=ρ on ∂(Ω^1,2_m ×Q³_m)×Ω³, t >0, (5.8) σ=ρinit in ΩM ×Q³_m, t= 0. (5.9) Boundary condition (5.8) reflects continuity of pressure between the blocks and the fractures. Let λ1 and λ2 be the first and second rows of the matrixλ. Then fm,f_M¹, and f_M² are given by

f_m(x, t) =− 1

|Q³| Z

Q³_m

φ(x₁, x₂, y₃)σ_t(x, y₃, t)dy₃, f_M¹(x, t) = 1

|Q³| Z

Q³_m

λ₁· ∇x1,x2,y3σ(x, y₃, t)dy₃,

f_M²(x, t) = 1

|Q³| Z

Q³_m

λ₂· ∇_x1,x2,y3σ(x, y₃, t)dy₃.

Together, the f_m, f_M¹, and f_M² terms are a fluid source, representing the total average flow out of the matrix blocks into the fracture sheet. Thef_m term alone contributes too much fluid because it takes into account lateral flow in the blocks as well as vertical. Since in Ω_M fluid is exchanged only through the top and bottom of each block, thef_M¹ andf_M² terms cancel the extra flow.

Theorem 5.1. The macroscopic model has a unique solution, (ρ, σ). Moreover, the solution of the microscopic model,(ρ^ε, σ^ε), converges to (ρ, σ) in the following weak sense: χ^ε_fΦ^∗ρ^ε*ΦρinH¹(J;L²(Ω)),χ^ε_fΛ∇ρ^ε*Λ^#ξ inL²(J;L²(Ω)), and fσ^ε* σinL²(Ω³;H¹(Ω^1,2_m ×Q³_m×J)).

The proof will be completed via the subsequent lemmas.

Lemma 5.2. For a subsequence of theε’s, we have the following weak convergence:

χ^ε_fΦ^∗ρ^ε*Φρ inH¹(J;L²(Ω)), (5.10) χ^ε_fΦ^∗ρ^ε*Φρ in H¹(ΩF×J), (5.11) χ^ε_fΦ^∗ρ^ε*Φρ inL²(Ω³;H¹(Ω^1,2_m ×J)), (5.12) χ^ε_fΛ∇ρ^ε*Λ^#ξ inL²(J;L²(Ω)), (5.13) fσ^ε* σ inL²(Ω³;H¹(Ω^1,2_m ×Q³_m×J)), (5.14) ρe^ε* τ3 inL²(Ω³;H¹(Ω^1,2_f ×Q³×J)), (5.15)

ρe^ε* τ in L²(Ω³;H¹(Ω^1,2_m ×Q³_f×J)). (5.16) The proof of the above convergence follow immediately from Lemma 3.2 and Corollary 4.4.

We denote the components of ξ as ξ= (ξ₁, ξ₂, ξ₃). From Corollary 4.4, we see that τ3=τ3(x, t) only, because Q³ is connected. Similarly, if we letτ1 and τ2 be the restrictions of τ to Ω_M ×Q^3,1_f ×J and Ω_M ×Q^3,2_f ×J respectively, we have τ1=τ1(x, t) andτ2=τ2(x, t) as well. Next we show thatρe^ε actually converges to ρ.

Lemma 5.3. τ₃=ρin L²(Ω³;H¹(Ω^1,2_f ×J))and ^τ1+τ₂ ² =ρin L²(Ω³;H¹(Ω^1,2_m × J)).

Proof. Letϕ∈C₀^∞(Ω_F×J). Then using Lemmas 4.2 and 4.3 we have (ρe^ε, ϕ)Ω_F×Q³×J= (]χ^ε_fρ^ε, ϕ)Ω_F×Q³×J = (χ^ε_fρ^ε,ϕ)e Ω_F×Q³×J

→( Φ

Φ^∗ρ, ϕ)Ω_F×Q³×J =|Q³|(ρ, ϕ)ΩF×J.

We also have (ρe^ε, ϕ)Ω_F×Q³×J → |Q³|(τ3, ϕ)ΩF×J, so the first result is clear. Simi- larly, forϕ∈C₀^∞(Ω_M ×J), Lemmas 4.2 and 4.3 give

(ρe^ε, ϕ)_ΩM×Q3

f×J = (]χ^ε_fρ^ε, ϕ)_ΩM×Q3×J = (χ^ε_fρ^ε,ϕ)e _ΩM×Q3×J

→(Φ

Φ^∗ρ, ϕ)_ΩM×Q3×J =|Q³_f|(ρ, ϕ)ΩM×J. Then the second result follows because

(ρe^ε, ϕ)_ΩM×Q³

f×J →(τ, ϕ)_ΩM×Q³

f×J=|Q^3,1_f |(τ1, ϕ)Ω_M×J+|Q^3,2_f |(τ2, ϕ)Ω_M×J

=|Q³_f|(τ1+τ2

2 , ϕ)Ω_M×J.

We claim that, in fact,τ₁=τ₂ and thus,τ=ρinL²(Ω³;H¹(Ω^1,2_m ×J)). Recall that the standard cellQ was constructed with the block centered in the fracture, that is,|Q^3,1_f |=|Q^3,2_f |= ¹₂|Q³_f|. If instead, we define a new standard cell ˇQ, with corresponding subsets, and repeat all the analysis, replacingτ1andτ2 withγ1 and γ2, then we have the following: τ1 = γ1 in L²(Ω³;H¹(Ω^1,2_m ×J)) and τ2 =γ2 in L²(Ω³;H¹(Ω^1,2_m ×J)). In addition, calculations similar to the ones in the proof of Lemma 5.3 give ^|Qˇ

3,1

f |γ1+|Qˇ^3,2_f |γ2

|Qˇ³_f| =ρin L²(Ω²;H¹(Ω¹_m×J)). This implies that τ1=τ2

|Qˇ^3,2_f | − |Q^3,2_f |

|Qˇ^3,1_f | − |Q^3,1_f |

,

which gives τ1 =τ2, since ^|Qˇ

3,2 f |−|Q^3,2_f |

|Qˇ^3,1_f |−|Q^3,1_f | = 1. Now we show that the weak limits ρ andσare solutions to the macroscopic model.

Lemma 5.4. The weak limit σ is a solution of (5.7), the block equation of the macroscopic model.

Proof. First takeψ∈L²(Ω³;L²(J;H₀¹(Ω^1,2_m ×Q³_m))) and define ψ(x¯ ₁, x₂, x₃, z₃, t) =

(ψ x₁, x₂, x₃,^z3−c_ε^ε(x3), t

ifz₃∈εQ³_m+c^ε(x₃)

0 otherwise.

Then, we have ¯ψ∈L²(Ω³;L²(J;H₀¹(Ω^ε_m))). We use ¯ψ as the test function in the block equation of the weak form of the microscopic model (3.11). So, for almost every fixedx₃∈Ω³, we have

(φ^εσ_t^ε,ψ(x¯ 3))Ω^ε_m×J+ (bλ^ε∇x₁,x₂,z₃σ^ε,∇x₁,x₂,z₃ψ(x¯ 3))Ω^ε_m×J = (f,ψ(x¯ 3))Ω^ε_m×J. Since ¯ψ= 0 in all blocks except the one that xis in, this simplifies to

(φ^εσ_t^ε,ψ(x¯ 3))_Ω^1,2

m ×(εQ³_m+c^ε(x₃))×J

+ (bλ^ε∇x₁,x₂,z₃σ^ε,∇x₁,x₂,z₃ψ(x¯ 3))_Ω1,2

m ×(εQ³_m+c^ε(x3))×J

= (f,ψ(x¯ ₃))_Ω1,2

m×(εQ³_m+c^ε(x₃))×J.

Next, we lety3= ^z3−c_ε^ε(x3) and use theQ³-periodicity ofφandbλto obtain (φfσ_t^ε, ψ)_Ω^1,2

m ×Q³_m×J+ (bλ∇^1/ε_x1,x2,y3σf^ε,∇^1/ε_x1,x2,y3ψ)_Ω^1,2

m×Q³_m×J= (f , ψ)e _Ω^1,2

m×Q³_m×J, where∇^1/εx1,x2,y3 = (_∂x^∂

1,_∂x^∂

2,¹_ε∂y^∂

3). Simplifying and integrating inx3, we get (φfσ_t^ε, ψ)Ω_M×Q³_m×J+ (λ∇x₁,x₂,y₃fσ^ε,∇x₁,x₂,y₃ψ)Ω_M×Q³_m×J = (f , ψ)e Ω_M×Q³_m×J. Finally, lettingε→0 gives us

(φσ_t, ψ)_ΩM×Q3

m×J+ (λ∇_x1,x2,y3σ,∇_x1,x2,y3ψ)_ΩM×Q3

m×J = (f, ψ)_ΩM×Q3

m×J. (5.17)

This is a weak form of (5.7).

Lemma 5.5. The weak limitρis a solution of (5.1)-(5.3),(5.5)-(5.6), the fracture equations of the macroscopic model.

Proof. Starting with (3.10), we let ε → 0 and consider the convergence of each term. The first and second terms are straightforward:

(Φ^∗ρ^ε_t, ϕ)Ω^ε_f×J= (χ^ε_fΦ^∗ρ^ε_t, ϕ)_Ω×J→(Φρt, ϕ)_Ω×J, (Λ∇ρ^ε,∇ϕ)Ω^ε_f×J = (χ^ε_fΛ∇ρ^ε,∇ϕ)_Ω×J →(Λ^#ξ,∇ϕ)_Ω×J.

Using the periodicity ofφand Lemmas 4.2-4.3, the third term converges as follows:

(φ^εσ_t^ε, ϕ)Ω^ε_m×J = 1

|Q³|(φfσ^ε_t,ϕ)e _ΩM×Q3 m×J

→ 1

|Q³|(φσt, ϕ)Ω_M×Q³_m×J =−(fm, ϕ)Ω_M×J. Lettingλ^ε_i,i= 1,2,3, be the rows ofλ^ε, we expand the fourth term,

(bλ^ε∇σ^ε,∇ϕ)Ω^ε_m×J

= (λ^ε₁· ∇^εσ^ε, ∂

∂x₁ϕ)Ω^ε_m×J+ (λ^ε₂· ∇^εσ^ε, ∂

∂x₂ϕ)Ω^ε_m×J+ε(λ^ε₃· ∇^εσ^ε, ∂

∂x₃ϕ)Ω^ε_m×J. The last term has an extraεand thus converges to 0 by Lemma 3.2. The conver- gence of other terms is shown using Lemmas 4.1-4.3:

(λ^ε_i∇^εσ^ε, ∂

∂xi

ϕ)_Ωε

m×J= 1

|Q³|(fλ^ε_i∇]^εσ^ε,^∂

∂xi

ϕ)_ΩM×Q3 m×J