\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2010(2010), No. 94, pp. 1--19.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2010 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2010/94\hfil A model for single-phase flow]
{A model for single-phase flow in layered porous media}

\author[D. J. Coffield Jr., A. M. Spagnuolo \hfil EJDE-2010/94\hfilneg]
{Daniel J. Coffield Jr., Anna Maria Spagnuolo}  % in alphabetical order

\address{Daniel J. Coffield Jr. \newline
University of Michigan-Flint\\
Mathematics Department\\
Flint, MI 48502-1950, USA}
\email{dcoffiel@umflint.edu, tel: (810) 762-3005, fax: (810) 766-6880}

\address{Anna Maria Spagnuolo \newline
Oakland University\\
Department of Mathematics and Statistics\\
Rochester, MI 48309-4485, USA}
\email{spagnuol@oakland.edu}

\thanks{Submitted February 23, 2010. Published July 13, 2010.}
\subjclass[2000]{76S05, 35B27}
\keywords{Layered media; naturally-fractured media; \hfill\break\indent
double-porosity model; homogenization}

\begin{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 $\varepsilon$, based on the
 fracture direction in the reservoir.  The parameter
 $\varepsilon$ 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.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}\label{sec1}

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 associated 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
\cite{BZK60,K69,Pirson,WR63} using physical arguments. Since
then, a more general form of a double-porosity model for single-phase
flow in totally fractured reservoirs was derived on physical grounds and
by using formal homogenization in
\cite{ArDoHu91,Ar89,da90,DPLAS87}, and then proven rigorously using
homogenization techniques \cite{ArDH90}. For the main ideas of fluid flow 
in porous media and
homogenization theory, see \cite{Be72,SaPa}, respectively.
In addition, a general treatment of homogenization and porous
media can be found in \cite{HoBook}.

Double-porosity models have also been used to better model partially
fissured media \cite{DoPeSh}. Additionally, the dual-porosity concept has been
applied to layered porous media. In \cite{DoAbPLHeNu} 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 \cite{ClSh} for a
porous medium with relatively thin rock cells. Furthermore,
multiple-porosity models have been derived using homogenization
techniques in \cite{ArDH90,DKPLS98,ShSpWr,SpW1,SpW2}.

Fractured reservoirs are often classified according to the extent and
characteristics 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 necessarily
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 \cite{ArDH90}, 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 \cite{DoAbPLHeNu}. 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 \ref{sec2}. In Section \ref{sec3}, we describe the
microscopic model. Section \ref{sec4} includes many technical lemmas and
Section \ref{sec5} defines the macroscopic limit model and includes the
details of how the microscopic model converges to the limit model.
Finally, Section \ref{sec6} briefly discusses some extensions to
this work.


\section{Preliminaries}\label{sec2}

We consider the reservoir
$\Omega=\Omega^1\times\Omega^2\times\Omega^3\subset \mathbb{R}^3$
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 by $Q$ and $\Omega$ is the interior of
the union of translations of $\overline{Q}$:
\[
\overline{\Omega}=\cup_{\alpha\in I(1)}(\overline{Q}+c_{\alpha}),
\]
where the $c_{\alpha}$'s are translations and $I(1)=\{0,1,2,\dots
,N(1)\}$ is an index set (See Fig. \ref{fig1}).

\begin{figure}[ht]
  \begin{center}
  \includegraphics[scale=.35]{fig1.jpg}
  \end{center}
\caption{The reservoir $\Omega=\Omega^1\times\Omega^2\times\Omega^3$.}
\label{fig1}
\end{figure}

Additionally, we require the cells to be disjoint:
$(Q+c_{\alpha})\cap(Q+c_{\beta})=\emptyset$ for $\alpha\neq\beta$.
Let $C(1)=\{c_{\alpha}:\alpha\in I(1)\}$ be the set of
translations. Then $(0,0,0)\in C(1)$ and we call it $c_0$.  Also,
there is a $c_{\alpha}=(0,0,r)\in C(1)$, call it $c_1$, with
$r\in\mathbb{R}$, such that $C(1)=\{\alpha c_1: \alpha\in I(1)\}$.
That is, $C(1)$ is made up of integer multiples of the standard
translation, $c_1$, and the translations only occur in the third
component.

The standard rectangular cell $Q=Q^1\times Q^2\times Q^3$ consists
of three parts: the matrix block, denoted by $Q_m$; the fracture
surrounding the matrix, denoted by $Q_f$; and the internal
boundary between $Q_m$ and $Q_f$, denoted by $\partial Q_m$ (See
Fig. \ref{fig2}).   The domain $Q_m$ is itself rectangular,
$Q_m=Q_m^1\times Q_m^2\times Q^3_m$, and is compactly contained in
$Q$.  The domain $Q_f$ is connected and $\partial Q_m$ is
Lipschitz and piecewise smooth.  We also define
$Q_f^1=Q^1\backslash \overline{Q_m^1}$, $Q_f^2=Q^2\backslash
\overline{Q_m^2}$,  and $Q_f^3=Q^3\backslash \overline{Q_m^3}$.
The domain $Q_f^3$ consists of two pieces, which we denote by
$Q_f^{3,1}$ and $Q_f^{3,2}$, with
$|Q_f^{3,1}|=|Q_f^{3,2}|=\frac{1}{2}|Q^3_f|$.

\begin{figure}[ht]
  \begin{center}
  \includegraphics[scale=.5]{fig2.jpg}
  \end{center}
\caption{The standard cell $Q$ and its subsets.}\label{fig2}
\end{figure}

To homogenize, we will shrink the cells in the direction of
periodicity, the vertical direction.  We introduce the notation,
${}^{\varepsilon}\!Q=Q^1\times Q^2\times \varepsilon  Q^3$, and
call this the $\varepsilon$-cell.  Similarly, we define
${}^{\varepsilon}\!Q_m = Q_m^1\times Q_m^2\times\varepsilon Q_m^3$
and ${}^{\varepsilon}\!Q_f = \{(y_1,y_2, \varepsilon y_3):
(y_1,y_2,y_3) \in Q_f\}$.  When $\varepsilon=1$, we have
${}^{\varepsilon}\!Q=Q$, the standard cell.  As $\varepsilon$
tends to $0$ the cells shrink linearly in the vertical direction
only.  For each $\varepsilon$, we have a new index set
$I(\varepsilon)=\{0,1,2,\dots , N(\varepsilon)\}$ and set of
translations $C(\varepsilon)=\{c_{\alpha}:\alpha\in
I(\varepsilon)\}$ so that the $\varepsilon$-cells fill up $\Omega$
(See Fig. \ref{fig3}).

\begin{figure}[ht]
  \begin{center}
 \includegraphics[scale=.40]{fig3.jpg}
  \end{center}
\caption{Vertical cross-sections through the middle of $\Omega$ showing 
the vertical homogenization as $\varepsilon\to 0$ and the 
$\varepsilon$-dependent subsets of $\Omega$.}\label{fig3}
\end{figure}

Next, we define a matrix domain and fracture domain for each
$\varepsilon >0$ as follows:
\[
\Omega^{\varepsilon}_m=\Omega\cap\cup_{\alpha\in I(\varepsilon)}
\left({}^{\varepsilon}\!Q_m + \varepsilon c_{\alpha}\right) \quad
\text{and} \quad \Omega^{\varepsilon}_f=\Omega\cap\cup_{\alpha\in
I(\varepsilon)} \left({}^{\varepsilon}\!Q_f + \varepsilon
c_{\alpha}\right) ,
\]
(See Fig. \ref{fig3}).  Note that the translations are $0$ in the
first two components so the $\varepsilon$ only affects the third
component of the translation $c_{\alpha}$. For convenience we
assume that the sequence of $\varepsilon$'s are such that
$\partial\Omega\subset\partial\Omega^{\varepsilon}_f$.

Next, we introduce notation for various components of $\Omega$.
Since, $\Omega^1=Q^1$ and $\Omega^2=Q^2$, we similarly define the
following: $\Omega^1_m=Q_m^1$, $\Omega^1_f=Q_f^1$, $\Omega^2_m=Q_m^2$,
and $\Omega^2_f=Q_f^2$.  For simplicity, we define
$\Omega^{1,2}=\Omega^1\times\Omega^2$,
$\Omega^{1,2}_m=\Omega^1_m\times\Omega^2_m$, and
$\Omega^{1,2}_f=(\Omega^1_f\times\Omega^2)\cup(\Omega^1\times\Omega^2_f)$.
 To describe the setting of the macroscopic model, we decompose
$\Omega$ into three pieces: the fracture sheet, denoted by
$\Omega_M=\Omega^{1,2}_m\times\Omega^3$, corresponding to
the part of $\Omega$ with matrix blocks in it; the side fractures,
denoted by $\Omega_F=\Omega^{1,2}_f\times\Omega^3$, that do not
shrink as $\varepsilon$ 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
$\partial\Omega^{1,2}_m\times\Omega^3$.  Note that
$\Omega^{1,2}_m\times \Omega_m^{3,\varepsilon}=\Omega^{\varepsilon}_m$
and $\Omega^{1,2}_m\times\Omega_f^{3,\varepsilon}
= \Omega^{\varepsilon}_f\backslash\overline{\Omega_F}$,
where $\Omega_m^{3,\varepsilon}$ and $\Omega_f^{3,\varepsilon}$
are two $\varepsilon$-dependent subsets of $\Omega^3$
(See Fig. \ref{fig4}).

\begin{figure}[ht]
  \begin{center}
  \includegraphics[scale=.40]{fig4.jpg}
  \end{center}
\caption{A vertical cross-section through the middle of $\Omega$ showing 
the $\varepsilon$-dependent subsets of $\Omega^3$.}\label{fig4}
\end{figure}

As the cells shrink, their horizontal components remain unchanged.
That is, $Q^1$, $Q^2$, and their subsets are not shrinking at all.
So, while the vertical components of $\Omega$ and $Q$ are on different
scales as $\varepsilon$ tends to $0$, the horizontal components
remain on the same scale and are, in fact,
the same: $Q^i=\Omega^i$, $Q_m^i=\Omega_m^i$, and $Q_f^i=\Omega_f^i$,
for $i=1,2$.  Because of this, we ignore the redundant spaces for
the horizontal components and refer only to $\Omega^1$, $\Omega^2$,
and their subsets.

For each $0<\varepsilon\leq 1$, we define a
function $c^{\varepsilon}:\Omega^3\to C(\varepsilon)$
by letting $c^{\varepsilon}(x_3)$ be the translation corresponding
to the $\varepsilon$-cell that contains $\Omega^{1,2}\times\{x_3\}$.
The function $c^{\varepsilon}$ can be properly defined by assigning
the $\varepsilon$-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 use $J$ to denote the open time
interval $(0,T)$.


\section{The Microscopic Model}\label{sec3}

In this section we pose the model on the microscopic scale.  For
$x=(x_1,x_2,x_3)$, let $\rho^{\varepsilon}(x,t)$ and
$\sigma^{\varepsilon}(x,t)$ be the fluid density in
$\Omega^{\varepsilon}_f$ and $\Omega^{\varepsilon}_m$
respectively.  Next, we assume the fluid is a liquid with
Newtonian viscosity, $\mu$, and that it satisfies the following
equations of state: $d\rho^\varepsilon= c\rho^\varepsilon
dp^\varepsilon$ and $d\sigma^\varepsilon= c\sigma^\varepsilon
dp^\varepsilon$, where $c$ is the constant (small) compressibility
and $p^\varepsilon$ is the pressure.

Let $\phi(x)$ and $k(x)$ be the porosity and permeability matrix,
respectively, in the blocks of the original reservoir,
$\Omega_m^{\varepsilon=1}$.  Since the blocks are periodically
distributed over $\Omega$ in the vertical direction, we assume
that $\phi$ and $k$ are vertically periodic with a period of $
Q^3$.  We then let
$\phi^{\varepsilon}(x)=\phi(x_1,x_2,\frac{x_3}{\varepsilon})$ and
$k^{\varepsilon}(x)=k(x_1,x_2,\frac{x_3}{\varepsilon})$ be the
porosity and permeability matrix in $\Omega^{\varepsilon}_m$ so
that $\phi^{\varepsilon}$ and $k^{\varepsilon}$ are $\varepsilon
Q^3$-periodic in the vertical direction.

Following \cite{ArDH90}, the fluid flow is assumed to be  governed
by Darcy's law in $\Omega^{\varepsilon}_m$ and
$\Omega^{\varepsilon}_f$.  We define $\Phi^*(x)$ and $K^*(x)$ to
be the porosity and scalar permeability of the original fracture
domain, defining them first over all of $\Omega$ and then
considering their restrictions to $\Omega_f^{\varepsilon=1}$. Now,
$\phi$, $\Phi^*$, $k$ and $K^*$ are all smooth and bounded, with
$\phi$, $\Phi^*$, and $K^*$ being uniformly positive.    Also, $k$
is 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 matrix $k^{\varepsilon}$. We
use a scaling similar to the one done in \cite{DoAbPLHeNu}. Let
$H$ be the  matrix
\[
H = \begin{pmatrix}
1 & 0 & 0  \\
0 & 1 & 0 \\
0 & 0 & \varepsilon
\end{pmatrix}.
\]
We then define a scaled permeability matrix,
$\widehat{k}^{\varepsilon}$, by
$\widehat{k}^{\varepsilon}= Hk^{\varepsilon}H$.  So,
\begin{align*}
\widehat{k}^{\varepsilon}(x)
&=  \begin{pmatrix}
k^{\varepsilon}_{11}(x) & k^{\varepsilon}_{12}(x)
 & \varepsilon k^{\varepsilon}_{13}(x)  \\
k^{\varepsilon}_{21}(x) & k^{\varepsilon}_{22}(x)
 & \varepsilon k^{\varepsilon}_{23}(x) \\
\varepsilon k^{\varepsilon}_{31}(x)
 & \varepsilon k^{\varepsilon}_{32}(x)
 & \varepsilon^2 k^{\varepsilon}_{33}(x)
\end{pmatrix} \\
& =  \begin{pmatrix}
k_{11}(x_1,x_2,\frac{x_3}{\varepsilon})
 & k_{12}(x_1,x_2,\frac{x_3}{\varepsilon})
 & \varepsilon k_{13}(x_1,x_2,\frac{x_3}{\varepsilon})  \\
k_{21}(x_1,x_2,\frac{x_3}{\varepsilon})
 & k_{22}(x_1,x_2,\frac{x_3}{\varepsilon})
 & \varepsilon k_{23}(x_1,x_2,\frac{x_3}{\varepsilon}) \\
\varepsilon k_{31}(x_1,x_2,\frac{x_3}{\varepsilon})
 & \varepsilon k_{32}(x_1,x_2,\frac{x_3}{\varepsilon})
 & \varepsilon^2 k_{33}(x_1,x_2,\frac{x_3}{\varepsilon})
\end{pmatrix} .
\end{align*}

Throughout we use $\nu$ 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
$\Omega^{\varepsilon}_f$ and $\Omega^{\varepsilon}_m$ can be
described in the standard way as follows:
\begin{gather}
\Phi^*\rho_t^{\varepsilon} - \nabla\cdot\left(\frac{K^*}{\mu c}
 \nabla\rho^{\varepsilon}\right)
 = f \quad  \text{in }\Omega^{\varepsilon}_f,\ t>0,\label{eq3.1}
\\
\frac{K^*}{\mu c}\nabla\rho^{\varepsilon}\cdot \nu
= \frac{\widehat{k}^{\varepsilon}}{\mu
c}\nabla\sigma^{\varepsilon}\cdot \nu \quad \text{on }
\partial \Omega^{\varepsilon}_m, \; t>0, \label{eq3.2}
\\
\frac{K^*}{\mu c}\nabla \rho^{\varepsilon}\cdot \nu = 0 \quad
 \text{on }  \partial\Omega,\; t > 0, \label{eq3.4}
\\
\rho^{\varepsilon} = \rho_{\rm init} \quad  \text{in }
 \Omega^{\varepsilon}_f,\; t = 0, \label{eq3.3}
\\
\phi^{\varepsilon}\sigma_t^{\varepsilon}
 - \nabla\cdot\left(\frac{\widehat{k}^{\varepsilon}}{\mu c}
 \nabla\sigma^{\varepsilon}\right) = f \quad
 \text{in }\Omega^{\varepsilon}_m,\; t>0, \label{eq3.5}
\\
\sigma^{\varepsilon} = \rho^{\varepsilon} \quad  \text{on }
  \partial \Omega^{\varepsilon}_m,\; t>0, \label{eq3.6}
\\
\sigma^{\varepsilon} = \rho_{\rm init} \quad  \text{in }
\Omega^{\varepsilon}_m,\ t = 0. \label{eq3.7}
\end{gather}
The function $f=f(x,t)$ represents external sources and sinks. The
boundary conditions \eqref{eq3.2} and \eqref{eq3.6} respectively
represent conservation of mass and continuity of pressure between
$\Omega^{\varepsilon}_f$ and $\Omega^{\varepsilon}_m$. Also,
\eqref{eq3.4} gives no-flow across the external boundary
$\partial\Omega$.  For the sake of convenience, we are ignoring
gravity in the model but could easily account for it.

The need for the $\varepsilon$-scaling in
$\widehat{k}^{\varepsilon}$ can be seen by considering
matrix-to-fracture flow as the blocks shrink. The scaling causes
the matrix blocks to be less permeable as
$\varepsilon\to0$, 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 \cite{ArDoHu91}, we choose the
$\varepsilon$-scaling in $\widehat{k}^{\varepsilon}$ to conserve
matrix-to-fracture flow in a certain sense. The original
($\varepsilon=1$) matrix-to-fracture flow is given by:
 \begin{equation} \label{eq3.9}
\begin{aligned}
\int_{\partial\Omega_m^{\varepsilon=1}} \frac{k(x)}{\mu c}
\nabla\sigma^{\varepsilon =1}(x)\cdot\nu\,ds(x)
& = \int_{\Omega_m^{\varepsilon=1}} \nabla\cdot\left(\frac{k(x)}{\mu
c} \nabla\sigma^{\varepsilon =1}(x)\right) \,dx \\
&= \sum_{i=0}^{N(1)}\int_{Q_m+c_i} \nabla\cdot
\left(\frac{k(x)}{\mu c} \nabla\sigma^{\varepsilon =1}(x)\right) dx ,
\end{aligned}
\end{equation}
where $N(1)$ is the number of cells in the original ($\varepsilon=1$)
domain.  Now, for arbitrary $\varepsilon<1$, we consider the
matrix-to-fracture flow without scaling.  Making a change of
variables and recalling that $k$ is $ Q^3$-periodic, the flow
is given by
\begin{equation} \label{eq3.10}
\begin{aligned}
&\int_{\partial\Omega^{\varepsilon}_m}
\frac{k^{\varepsilon}(x)}{\mu c}
\nabla\sigma^{\varepsilon}(x)\cdot\nu\,ds(x)\\
&= \int_{\Omega^{\varepsilon}_m} \nabla\cdot\left(\frac{k^{\varepsilon}(x)}{\mu c} 
\nabla\sigma^{\varepsilon}(x)\right) dx \\
&= \sum_{i=0}^{N(\varepsilon)}\quad \int_{{}^{\varepsilon}\!Q_m+\varepsilon c_i} 
\nabla\cdot\left(\frac{k^{\varepsilon}(x)}{\mu c} \nabla\sigma^{\varepsilon}(x)\right) dx \\
& = \varepsilon \sum_{i=0}^{N(\varepsilon)}\quad \int_{Q_m+ c_i}
\nabla^{1/\varepsilon}_z\cdot\left(\frac{k(z_1,z_2,z_3)}{\mu c}
\nabla^{1/\varepsilon}_z\sigma^{\varepsilon}(z_1,z_2,\varepsilon
z_3)\right) dz,
\end{aligned}
\end{equation}
where $\nabla^{1/\varepsilon}_z$ means
$\big(\frac{\partial}{\partial z_1},\frac{\partial}{\partial z_2},
\frac{1}{\varepsilon}\frac{\partial}{\partial z_3}\big)$.
Comparing \eqref{eq3.9} and \eqref{eq3.10}, we see that the
extra $\varepsilon$ is accounted for since $N(\varepsilon)$
is on the order of $\varepsilon^{-1}N(1)$.  But, the extra
$\frac{1}{\varepsilon}$'s are not accounted for.
By scaling $k^{\varepsilon}$ as we have done in
$\widehat{k}^{\varepsilon}$, we account for them and keep the
flow from blowing up.  Thus, choosing
$\widehat{k}^{\varepsilon}=Hk^{\varepsilon}H$ seems to be an
appropriate choice for conserving the matrix-to-fracture flow
as we homogenize.

To simplify things, we introduce more notation:
\begin{gather*}
\Phi(x) = \begin{cases}
\Phi^*(x) & \text{if }  x \in \Omega_F\\
\frac{|Q^3_f|}{|Q^3|}\Phi^*(x) & \text{if } x\in \Omega_M ,
\end{cases}
\quad
\Lambda (x) = \Lambda_F (x) = \frac{K^*(x)}{\mu c},\\
\Lambda_M (x) = \frac{|Q^3_f|}{|Q^3|}\frac{K_M(x)}{\mu c},\quad
\lambda^{\varepsilon} (x) = \frac{k^{\varepsilon}(x)}{\mu c}, \quad
\widehat{\lambda}^{\varepsilon} (x) =
\frac{\widehat{k}^{\varepsilon}(x)}{\mu c}.
\end{gather*}
Throughout we assume that $f\in L^2(J;L^2(\Omega))$,
$\rho_{\rm init}\in L^2(\Omega)$, and use $(\cdot,\cdot)_X$ to
denote the standard inner product of $L^2(X)$.
Then for $\varphi \in L^2(J;H^1(\Omega))$, we use the boundary
conditions and \eqref{eq3.5} to see that
\begin{align*}
- (\nabla\cdot(\Lambda\nabla\rho^{\varepsilon}),\varphi)_{\Omega^{\varepsilon}_f}
&= (\Lambda\nabla\rho^{\varepsilon},\nabla\varphi)_{\Omega^{\varepsilon}_f} 
- (\Lambda\nabla\rho^{\varepsilon}\cdot\nu_f,\varphi)_{\partial\Omega^{\varepsilon}_m} 
- (\Lambda\nabla\rho^{\varepsilon}\cdot\nu,\varphi)_{\partial\Omega}\\
&= (\Lambda\nabla\rho^{\varepsilon},\nabla\varphi)_{\Omega^{\varepsilon}_f} 
+(\widehat{\lambda}^{\varepsilon}\nabla\sigma^{\varepsilon}\cdot\nu_m,
\varphi)_{\partial\Omega^{\varepsilon}_m}\\
&=(\Lambda\nabla\rho^{\varepsilon},\nabla\varphi)_{\Omega^{\varepsilon}_f} 
+(\phi^{\varepsilon}\sigma^{\varepsilon}_t,\varphi)_{\Omega^{\varepsilon}_m}
 -(f,\varphi)_{\Omega^{\varepsilon}_m}  + (\widehat{\lambda}^{\varepsilon}
 \nabla\sigma^{\varepsilon},\nabla\varphi)_{\Omega^{\varepsilon}_m},
\end{align*}
where $\nu_f$ and $\nu_m$ are the outward pointing unit
normals from $\Omega^{\varepsilon}_f$ and $\Omega^{\varepsilon}_m$
respectively.  So a weak form for the microscopic model is
\begin{gather}
\begin{aligned}
&(\Phi^*\rho_t^{\varepsilon}, \varphi)_{\Omega^{\varepsilon}_f\times J}
+ (\Lambda\nabla\rho^{\varepsilon},
\nabla\varphi)_{\Omega^{\varepsilon}_f\times J}
+ (\phi^{\varepsilon}\sigma^{\varepsilon}_t,
 \varphi)_{\Omega^{\varepsilon}_m\times J}
 + (\widehat{\lambda}^{\varepsilon}\nabla\sigma^{\varepsilon},
 \nabla\varphi)_{\Omega^{\varepsilon}_m\times J}\\
&= (f,\varphi)_{\Omega\times J}, \quad \varphi\in L^2(J;H^1(\Omega)) ,
\end{aligned} \label{eq5.1}
\\
(\phi^{\varepsilon}\sigma^{\varepsilon}_t,\psi)_{\Omega^{\varepsilon}_m\times J}
 +(\widehat{\lambda}^{\varepsilon}\nabla\sigma^{\varepsilon},
 \nabla\psi)_{\Omega^{\varepsilon}_m\times J}
= (f,\psi)_{\Omega^{\varepsilon}_m\times J},
\quad \psi\in L^2(J;H^1_0(\Omega^{\varepsilon}_m)),\label{eq5.2}
\\
\rho^{\varepsilon}=\sigma^{\varepsilon} \quad\text{for }
 x\in\partial\Omega^{\varepsilon}_m,\; t>0,\label{eq5.3}
\\
\rho^{\varepsilon}=\rho_{\rm init},\quad
\sigma^{\varepsilon}=\rho_{\rm init} \quad\text{for }
  x\in\Omega,\; t=0.\label{eq5.4}
\end{gather}

\begin{theorem}\label{thm1}
For each $\varepsilon$, there exists a unique solution,
$\rho^\varepsilon$ and $\sigma^\varepsilon$, to the system
\eqref{eq5.1}--\eqref{eq5.4}.  Additionally,
$\rho^\varepsilon\in H^1(J;
L^2(\Omega^{\varepsilon}_f))\cap L^2(J;H^1(\Omega^{\varepsilon}_f))$
 and $\sigma^\varepsilon\in H^1(J;L^2(\Omega^{\varepsilon}_m))\cap L^2
(J;H^1(\Omega^{\varepsilon}_m))$.
\end{theorem}

\begin{proof}
First, let $\chi_f^{\varepsilon}$ be the characteristic function
of $\Omega^{\varepsilon}_f$ and
\begin{gather*}
\theta^{\varepsilon}(x) =
\begin{cases}
\rho^{\varepsilon} & \text{for }  x \in \Omega^{\varepsilon}_f\\
\sigma^{\varepsilon} & \text{for }  x \in \Omega^{\varepsilon}_m ,
\end{cases}
 \\
\alpha^{\varepsilon}=\chi_f^{\varepsilon}\Phi^*+
(1-\chi_f^{\varepsilon})\phi^{\varepsilon},\quad
\kappa^{\varepsilon}=\chi_f^{\varepsilon}\Lambda
+(1-\chi_f^{\varepsilon})\widehat{\lambda}^{\varepsilon}.
\end{gather*}
Then \eqref{eq5.1} and \eqref{eq5.4} can be rewritten as
\begin{align}
(\alpha^{\varepsilon}\theta_t^{\varepsilon}, \varphi)_{\Omega\times J} 
+ (\kappa^{\varepsilon}\nabla\theta^{\varepsilon}, \nabla\varphi)_{\Omega\times J} 
&= (f,\varphi)_{\Omega\times J}, \quad \forall \ \varphi\in L^2(J;H^1(\Omega)), \\
\theta^{\varepsilon} &= \rho_{\rm init} \quad\text{for}\ x\in\Omega,\
t=0  .
\end{align}
This is a single parabolic problem with discontinuous coefficients.
It is well known that this problem has a unique solution in
$H^1(J;L^2(\Omega))\cap L^2(J;H^1(\Omega))$.  By restricting
the solution to $\Omega^{\varepsilon}_f$ and $\Omega^{\varepsilon}_m$,
we obtain $\rho^\varepsilon$ and $\sigma^\varepsilon$ respectively.
\end{proof}

Hereafter we use $C$ to denote a positive constant that is
independent of $\varepsilon$, though not necessarily the same
one at each occurrence.  Also, we use $\nabla^{\varepsilon}$ to
denote $(\frac{\partial}{\partial x_1},\frac{\partial}{\partial x_2},
\varepsilon \frac{\partial}{\partial x_3})$.
Then, from the standard energy estimates used to prove
Theorem \ref{thm1}, we obtain the following result.

\begin{lemma}\label{lem1}
\begin{gather*}
\|\rho^{\varepsilon}\|_{L^2(J;L^2(\Omega^{\varepsilon}_f))}
 + \|\sigma^{\varepsilon}\|_{L^2(J;L^2(\Omega^{\varepsilon}_m))} \leq C, \\
\|\rho^{\varepsilon}_t\|_{L^2(J;L^2(\Omega^{\varepsilon}_f))}
 + \|\sigma^{\varepsilon}_t\|_{L^2(J;L^2(\Omega^{\varepsilon}_m))} \leq C, \\
\|\nabla\rho^{\varepsilon}\|_{L^2(J;L^2(\Omega^{\varepsilon}_f))}
+ \|\nabla^{\varepsilon}\sigma^{\varepsilon}\|_{L^2
(J;L^2(\Omega^{\varepsilon}_m))}
\leq C.
\end{gather*}
\end{lemma}

\section{Technical Lemmas} \label{sec4}

To prove the main result we will need to use various technical
lemmas which we state in this section.  First, for a measurable
function $\phi$ on $\Omega$, a subset of $\Omega$, or a subset
of the boundary of $\Omega$, we define
\[
\widetilde{\phi}(x_1,x_2,x_3,y_3) = \phi(x_1,x_2, \varepsilon y_3
+ c^{\varepsilon}(x_3)) .
\]
We call the operator mapping $\phi$ to $\widetilde{\phi}$
the dilation operator.  Recalling the definition of $c^{\varepsilon}$,
we are actually defining a family of operators, one for each
$\varepsilon$, $0<\varepsilon\leq 1$.  However, for convenience
we will only write ``\,$_{\widetilde{~}}$\,", and leave the
$\varepsilon$-dependence as implicitly understood.
We first note that the dilation operator commutes with addition
and multiplication:
$\widetilde{\phi\psi}= \widetilde{\phi}\widetilde{\psi}$ and
$\widetilde{\phi + \psi}= \widetilde{\phi} + \widetilde{\psi}$.
In the following lemmas, we prove some deeper properties of the
dilation operator that will be used extensively.

\begin{lemma}\label{lem7.1}
If $\psi \in H^1(\Omega)$, then
$\frac{\partial}{\partial y_3} \widetilde{\psi}
= \varepsilon \widetilde{\frac{\partial}{\partial x_3} \psi}$.
\end{lemma}

\begin{lemma}\label{lem7.3}
If $\psi, \varphi \in L^2(\Omega)$ (or $H^1(\Omega)$ when appropriate),
then for $r = m, f,$ or blank,
\begin{gather*}
(\widetilde{\psi},\widetilde{\varphi})_{\Omega^{1,2}_{p(r)}\times
\Omega^3\times Q^3_r} =
|Q^3|(\psi,\varphi)_{\Omega^{1,2}_{p(r)}\times \Omega^{3,\varepsilon}_r},\\
\|\widetilde{\psi}\|_{L^2(\Omega^{1,2}_{p(r)}\times
\Omega^3\times Q^3_r)} =
|Q^3|^{1/2}\|\psi\|_{L^2(\Omega^{1,2}_{p(r)}\times
\Omega^{3,\varepsilon}_r)},\\
\|\frac{\partial}{\partial y_3}
\widetilde{\psi}\|_{L^2(\Omega^{1,2}_{p(r)}\times
\Omega^3\times Q^3_r)} =
\varepsilon |Q^3|^{1/2}\|\frac{\partial}{\partial x_3}
\psi\|_{L^2(\Omega^{1,2}_{p(r)}\times \Omega^{3,\varepsilon}_r)}, \\
(\widetilde{\psi},\varphi)_{\Omega^{1,2}_r\times \Omega^3\times
Q^3} = (\psi,\widetilde{\varphi})_{\Omega^{1,2}_r\times
\Omega^3\times Q^3},
\end{gather*}
where $\Omega^{3,\varepsilon}$ means $\Omega^3$ and
$p(r) = \begin{cases}
m & \text{if }  r = f\text{ or } m,\\
f &  \text{if }  r = blank.
\end{cases}
$
\end{lemma}

\begin{proof}
Letting $z_3=\varepsilon y_3 + c^{\varepsilon}(x_3)$,
the first result is a calculation.
\begin{align*}
(\widetilde{\psi},\widetilde{\varphi})_{\Omega^{1,2}_{p(r)}\times 
\Omega^3\times Q^3_r}
& =  \int_{\Omega^{1,2}_{p(r)}} \int_{\Omega^3} \int_{Q^3_r} \psi(x_1,x_2,
\varepsilon y_3+c^\varepsilon (x_3)) \varphi(x_1,\varepsilon y_3+c^\varepsilon 
(x_3)) \, dy_3\, dx\\
& = \int_{\Omega^{1,2}_{p(r)}} \int_{\Omega^3} \int_{\varepsilon
Q^3_r+c^\varepsilon (x_3)} \psi(x_1,x_2,z_3) \varphi(x_1,x_2,z_3)
\varepsilon^{-1} \, dz_3\, dx \\
& = |Q^3| \int_{\Omega^{1,2}_{p(r)}} \int_{\Omega^{3,\varepsilon}_r}
\psi(x_1,x_2,z_3) \varphi(x_1,x_2,z_3) \, dz_3\, dx_1 \, dx_2  \\
& = |Q^3|(\psi,\varphi)_{\Omega^{1,2}_{p(r)} \times
\Omega^{3,\varepsilon}_r}.
\end{align*}
The next two results follow from the first result and the previous
lemma.  The last result is also a calculation.
\begin{align*}
&(\widetilde{\psi},\varphi)_{\Omega^{1,2}_r\times \Omega^3\times Q^3}\\
&= \int_{\Omega^{1,2}_r} \int_{\Omega^3} \int_{Q^3}
 \psi(x_1,x_2,\varepsilon y_3+c^\varepsilon (x_3)) \varphi(x_1,x_2,x_3)
 \, dy_3\, dx \\
&= \int_{\Omega^{1,2}_r} \int_{\Omega^3}
 \int_{\varepsilon Q^3+c^\varepsilon (x_3)} \psi(x_1,x_2,z_3)
 \varphi(x_1,x_2,x_3) \varepsilon^{-1} \, dz_3 \,dx \\
&= \varepsilon^{-1} \int_{\Omega^{1,2}_r} \int_{\Omega^3}
 \psi(x_1,x_2,z_3) \Big(\int_{\varepsilon Q^3+c^\varepsilon(z_3) }
 \varphi(x_1,x_2,x_3) \,dx_3\Big) dz_3\, dx_1 \, dx_2  \\
&= \varepsilon^{-1} \int_{\Omega^{1,2}_r} \int_{\Omega^3}
 \psi(x_1,x_2,z_3) \Big(\int_{Q^3} \varphi(x_1,x_2,\varepsilon y_3
 +c^{\varepsilon}(z_3)) \varepsilon \,dy_3\Big) dz_3\, dx_1 \, dx_2  \\
&= (\psi,\widetilde{\varphi})_{\Omega^{1,2}_r\times
\Omega^3\times Q^3}.
\end{align*}
\end{proof}

\begin{lemma}\label{lem7.5}
Let $s=F$ or $M$ and $r=m$, $f$, or $blank$.
If $\varphi\in L^2(\Omega_s)$ is considered as a function
in $L^2(\Omega_s\times Q^3_r)$ that is constant in $y_3$,
then $\widetilde{\varphi}\to\varphi$ strongly in
$L^2(\Omega_s\times Q^3_r)$ as $\varepsilon\to 0$.
\end{lemma}

\begin{proof}
Using the Dominated Convergence Theorem, the result is straightforward
for $\psi \in C^{\infty}_0(\Omega_s)$:
\begin{align*}
&\lim_{\varepsilon\to 0} \int_{\Omega_s\times Q_r^2}
(\widetilde{\psi}(x,y_2)-\psi(x))^2 \,dx\,dy_2\\
&=\int_{\Omega_s\times Q_r^2} (\lim_{\varepsilon\to
0}\psi(x_1,\varepsilon y_2 + c^{\varepsilon}(x_2))
-\psi(x))^2 \,dx\,dy_2=0.
\end{align*}
The full result then follows from Lemma \ref{lem7.3} and the
fact that functions in $C^{\infty}_0(\Omega_s)$ are dense
in $L^2(\Omega_s)$.
\end{proof}

\begin{corollary}\label{cor1}
We have the following estimates:
\begin{gather*}
\|\widetilde{\sigma^{\varepsilon}}\|_{L^2(\Omega^3;H^1(\Omega^{1,2}_m
\times Q^3_m\times J))}
\leq C,  \\
\|\widetilde{\rho^{\varepsilon}}\|_{L^2(\Omega^3;H^1(\Omega^{1,2}_f
\times Q^3\times J))}
\leq C,  \\
\|\widetilde{\rho^{\varepsilon}}\|_{L^2(\Omega^3;H^1(\Omega^{1,2}_m
\times Q^3_f\times J))}
\leq C,  \\
\|\frac{\partial}{\partial y_3}\widetilde{\rho^{\varepsilon}}\|_{L^2(\Omega_F;
L^2( Q^3;L^2(J)))}
\leq \varepsilon C,  \\
\|\frac{\partial}{\partial
y_3}\widetilde{\rho^{\varepsilon}}\|_{L^2(\Omega_M;L^2(
Q^3_f;L^2(J)))} \leq \varepsilon C.
\end{gather*}
\end{corollary}

The proofs for the above estimates are straighforward calculations
using Lemmas \ref{lem1} and \ref{lem7.3}.

\section{The Main Result}\label{sec5}

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

\begin{figure}[ht]
  \begin{center}
  \includegraphics[scale=.4]{fig5.jpg}
  \end{center}
\caption{The macroscopic model domains, $\Omega_F$ and $\Omega_M$, including 
a horizontal cross-section.}\label{fig13}
\end{figure}

Let $\rho^F(x,t)$ and $\rho^M(x,t)$ be the macroscopic,
averaged fluid density in side fractures and fracture sheet,
$\Omega_F$ and $\Omega_M$, respectively (See Fig. \ref{fig13}).
 Also, let
\[
\rho(x,t) = \begin{cases}
\rho^F(x,t) & \text{if }  x \in \Omega_F ,\\
\rho^M(x,t) & \text{if }  x\in \Omega_M ,
\end{cases}
\quad \Lambda^{\#}(x) = \begin{cases}
\Lambda_F(x) & \text{if }  x\in \Omega_F ,\\
\Lambda_M(x) & \text{if }  x\in \Omega_M .
\end{cases}
\]

We let $\nabla=(\frac{\partial}{\partial x_1},
\frac{\partial}{\partial x_2},\frac{\partial}{\partial x_3})$,
$\nabla_{x_1,x_2}=(\frac{\partial}{\partial x_1},
\frac{\partial}{\partial x_2})$, and
$\nabla_{x_1,x_2,y_3}=(\frac{\partial}{\partial x_1},
\frac{\partial}{\partial x_2},\frac{\partial}{\partial y_3})$.
Then, with $f_m$, $f^1_M$, and $f^2_M$ defined as below,
$\rho^F$ and $\rho^M$ satisfy the following:
\begin{gather}
\Phi\rho^F_t - \nabla\cdot\left(\Lambda_F\nabla\rho^F\right)
= f  \quad  \text{in }\Omega_F,\; t>0, \label{eq4.2}
\\
\Phi\rho^M_t - \nabla_{x_1,x_2}\cdot\left(\Lambda_M\nabla_{x_1,x_2}
\rho^M\right) - \frac{\partial}{\partial x_1}f^1_M
- \frac{\partial}{\partial x_2}f^2_M= f + f_m \quad
\text{in }\Omega_M,\; t>0, \label{eq4.3}
\\
\Lambda_F\nabla\rho^F\cdot\nu = 0 \quad  \text{on }
 \partial\Omega\backslash (\Omega^{1,2}_m\times\partial\Omega^3),\; t>0,
 \label{eq4.4}
\\
\rho = \rho_{\rm init} \quad  \text{in }\Omega,\; t = 0, \label{eq4.6}
\\
\frac{\partial}{\partial x_1}\rho^F
= \frac{|Q^3_f|}{|Q^3|}\frac{\partial}{\partial x_1}\rho^M, \quad
f^1_M=0 \quad  \text{on }
\partial\Omega^1_m\times\Omega^2_m\times\Omega^3,\; t>0, \label{eq4.7}
\\
\frac{\partial}{\partial x_2}\rho^F = \frac{|Q^3_f|}{|Q^3|}
\frac{\partial}{\partial x_2}\rho^M, \quad f^2_M =0 \quad  \text{on }
 \Omega^1_m\times\partial\Omega^2_m\times\Omega^3,\; t>0.
\label{eq4.7a}
\end{gather}

The boundary condition \eqref{eq4.4} holds on the boundary of
$\Omega$, excluding the top and bottom of $\Omega_M$.
The conditions \eqref{eq4.7} and \eqref{eq4.7a} hold on the
internal boundary between $\Omega_F$ and $\Omega_M$.
As $\varepsilon$ tends to zero, the blocks shrink to become
horizontal plates.  There are infinitely many of them, one for
each $x_3$.  Let $\sigma(x,y_3,t)$ be the macroscopic fluid density
in the matrix blocks and note that the macroscopic porosity and
scaled permeability matrix are $\phi(x_1,x_2,y_3)$ and
$\lambda(x_1,x_2,y_3)$ respectively.  Then $\sigma$ satisfies
the following:
\begin{gather}
\phi\sigma_t - \nabla_{x_1,x_2,y_3}\cdot
\left(\lambda\nabla_{x_1,x_2,y_3}\sigma\right) = f(x,t) \quad
\text{in }\Omega_M\times Q^3_m,\; t>0 \label{eq4.9},\\
\sigma = \rho \quad \text{on } \ \partial(\Omega^{1,2}_m\times
Q^3_m)\times\Omega^3,\; t>0 \label{eq4.10},\\
\sigma= \rho_{\rm init} \quad \text{in }\Omega_M\times Q^3_m,\; t = 0
\label{eq4.11}.
\end{gather}

Boundary condition \eqref{eq4.10} reflects continuity of pressure
between the blocks and the fractures. Let $\lambda_1$ and
$\lambda_2$ be the first and second rows of the matrix $\lambda$.
Then $f_m$, $f^1_M$, and $f^2_M$ are given by
\begin{gather*}
f_m(x,t) = -\frac{1}{|Q^3|}\int_{Q^3_m} \phi(x_1,x_2,y_3)
 \sigma_t(x,y_3,t)\,dy_3, \\
f^1_M(x,t) = \frac{1}{|Q^3|}\int_{Q^3_m} \lambda_1\cdot
\nabla_{x_1,x_2,y_3} \sigma(x,y_3,t)\,dy_3,\\
f^2_M(x,t) = \frac{1}{|Q^3|}\int_{Q^3_m} \lambda_2\cdot
\nabla_{x_1,x_2,y_3} \sigma(x,y_3,t)\,dy_3.
\end{gather*}
Together, the $f_m$, $f^1_M$, and $f^2_M$ terms are a fluid source,
representing the total average flow out of the matrix blocks
into the fracture sheet.  The $f_m$ term alone contributes too
much fluid because it takes into account lateral flow in the
blocks as well as vertical.  Since in $\Omega_M$ fluid is
exchanged only through the top and bottom of each block,
the $f^1_M$ and $f^2_M$ terms cancel the extra flow.

\begin{theorem}\label{thm2}
The macroscopic model has a unique solution, ($\rho,\sigma)$.
Moreover, the solution of the microscopic model,
$(\rho^\varepsilon,\sigma^\varepsilon)$, converges to
$(\rho,\sigma)$ in the following weak sense:
$\chi^{\varepsilon}_f \Phi^*\rho^{\varepsilon} \rightharpoonup
\Phi \rho$ in $H^1(J;L^2(\Omega))$,
$\chi^{\varepsilon}_f \Lambda\nabla\rho^{\varepsilon}
\rightharpoonup \Lambda^{\#} \xi$ in $L^2(J;L^2(\Omega))$,
and $\widetilde{\sigma^{\varepsilon}} \rightharpoonup \sigma$
in $L^2(\Omega^3;H^1(\Omega^{1,2}_m\times Q^3_m\times J))$.
\end{theorem}

The proof will be completed via the subsequent lemmas.

\begin{lemma}\label{lem3}
For a subsequence of the $\varepsilon$'s, we have the following
weak convergence:
\begin{gather}
\chi^{\varepsilon}_f \Phi^*\rho^{\varepsilon} \rightharpoonup
\Phi \rho \quad  \text{in } H^1(J;L^2(\Omega)), \label{eq9.1}\\
\chi^{\varepsilon}_f \Phi^*\rho^{\varepsilon} \rightharpoonup
\Phi \rho \quad  \text{in}\ H^1(\Omega_F\times J),  \label{eq9.4}\\
\chi^{\varepsilon}_f \Phi^*\rho^{\varepsilon} \rightharpoonup
\Phi \rho \quad  \text{in } L^2(\Omega^3;H^1(\Omega^{1,2}_m\times J)),
 \label{eq9.5}\\
\chi^{\varepsilon}_f \Lambda\nabla\rho^{\varepsilon} \rightharpoonup
\Lambda^{\#} \xi \quad   \text{in } L^2(J;L^2(\Omega)), \label{eq9.6}\\
\widetilde{\sigma^{\varepsilon}} \rightharpoonup \sigma\quad
 \text{in } L^2(\Omega^3;H^1(\Omega^{1,2}_m\times Q^3_m\times J)),
 \label{eq9.7}\\
\widetilde{\rho^{\varepsilon}} \rightharpoonup \tau_3 \quad \text{in }
 L^2(\Omega^3;H^1(\Omega^{1,2}_f\times Q^3\times J)),\label{eq9.8}\\
\widetilde{\rho^{\varepsilon}} \rightharpoonup \tau\quad
\text{in } L^2(\Omega^3;H^1(\Omega^{1,2}_m\times Q^3_f\times J)).
\label{eq9.9}
\end{gather}
\end{lemma}

The proof of the above convergence follow immediately from
Lemma \ref{lem1} and Corollary \ref{cor1}.


We denote the components of $\xi$ as $\xi=(\xi_1,\xi_2,\xi_3)$.
 From Corollary \ref{cor1}, we see that $\tau_3=\tau_3(x,t)$ only,
because $Q^3$ is connected.  Similarly, if we let $\tau_1$ and
$\tau_2$ be the restrictions of $\tau$ to
$\Omega_M\times Q^{3,1}_f \times J$ and
$\Omega_M\times Q^{3,2}_f \times J$ respectively, we have
$\tau_1=\tau_1(x,t)$ and $\tau_2=\tau_2(x,t)$ as well.
Next we show that $\widetilde{\rho^{\varepsilon}}$ actually
converges to $\rho$.

\begin{lemma}\label{lem9.3}
$\tau_3 = \rho$ in $L^2(\Omega^3;H^1(\Omega^{1,2}_f\times J))$ and
 $\frac{\tau_1+\tau_2}{2} = \rho$ in
$L^2(\Omega^3;H^1(\Omega^{1,2}_m\times J))$.
\end{lemma}

\begin{proof}
Let $\varphi \in C^{\infty}_0(\Omega_F\times J)$.  Then using
Lemmas \ref{lem7.3} and \ref{lem7.5} we have
\begin{align*}
(\widetilde{\rho^{\varepsilon}},\varphi)_{\Omega_F\times Q^3\times J}
&= (\widetilde{\chi^{\varepsilon}_f \rho^{\varepsilon}},
 \varphi)_{\Omega_F\times Q^3\times J}
= (\chi^{\varepsilon}_f \rho^{\varepsilon},\widetilde{\varphi}
 )_{\Omega_F\times Q^3\times J}\\
&\to (\frac{\Phi}{\Phi^*} \rho,\varphi)_{\Omega_F\times Q^3\times J}
=|Q^3|(\rho,\varphi)_{\Omega_F\times J}.
 \end{align*}
We also have $(\widetilde{\rho^{\varepsilon}},
 \varphi)_{\Omega_F\times Q^3\times J} \to |Q^3|(\tau_3,
 \varphi)_{\Omega_F\times J}$,
so the first result is clear.  Similarly, for
$\varphi \in C^{\infty}_0(\Omega_M\times J)$, Lemmas \ref{lem7.3}
and \ref{lem7.5} give
\begin{align*}
(\widetilde{\rho^{\varepsilon}},\varphi)_{\Omega_M\times
Q^3_f\times J}
&= (\widetilde{\chi^{\varepsilon}_f
\rho^{\varepsilon}},\varphi)_{\Omega_M\times Q^3\times J} =
(\chi^{\varepsilon}_f
\rho^{\varepsilon},\widetilde{\varphi})_{\Omega_M\times Q^3\times J}\\
&\to (\frac{\Phi}{\Phi^*}
\rho,\varphi)_{\Omega_M\times Q^3\times J} =
|Q^3_f|(\rho,\varphi)_{\Omega_M\times J}.
\end{align*}
Then the second result follows because
\begin{align*}
(\widetilde{\rho^{\varepsilon}},\varphi)_{\Omega_M\times
Q^3_f\times J} \to (\tau,\varphi)_{\Omega_M\times
Q^3_f\times J}& =  |Q^{3,1}_f|(\tau_1,\varphi)_{\Omega_M\times J} +
|Q^{3,2}_f|(\tau_2,\varphi)_{\Omega_M\times J}\\
 &= |Q^3_f|(\frac{\tau_1+\tau_2}{2},\varphi)_{\Omega_M\times J}.
\end{align*}
\end{proof}

We claim that, in fact, $\tau_1=\tau_2$ and thus, $\tau=\rho$
in $L^2(\Omega^3;H^1(\Omega^{1,2}_m\times J))$.
Recall that the standard cell $Q$ was constructed with the block
 centered in the fracture, that is,
$|Q_f^{3,1}|=|Q_f^{3,2}|=\frac{1}{2}|Q_f^3|$.
If instead, we define a new standard cell $\check{Q}$,
with corresponding subsets, and repeat all the analysis,
replacing $\tau_1$ and $\tau_2$ with $\gamma_1$ and $\gamma_2$,
then we have the following: $\tau_1=\gamma_1$ in
$L^2(\Omega^3;H^1(\Omega^{1,2}_m\times J))$ and $\tau_2=\gamma_2$
in $L^2(\Omega^3;H^1(\Omega^{1,2}_m\times J))$.  In addition,
calculations similar to the ones in the proof of Lemma \ref{lem9.3}
give $\frac{|\check{Q}^{3,1}_f|\gamma_1+|\check{Q}^{3,2}_f
|\gamma_2}{|\check{Q}^3_f|} = \rho$ in $L^2(\Omega^2;
H^1(\Omega^1_m\times J))$.  This implies that
\[
\tau_1=\tau_2\Big(\frac{|\check{Q}^{3,2}_f|-|Q^{3,2}_f|}{
|\check{Q}^{3,1}_f|-|Q^{3,1}_f|}\Big),
\]
which gives $\tau_1=\tau_2$, since
$\frac{|\check{Q}^{3,2}_f|-|Q^{3,2}_f|}{|\check{Q}^{3,1}_f|-|Q^{3,1}_f|}=1$.
Now we show that the weak limits $\rho$ and $\sigma$ are solutions
to the macroscopic model.


\begin{lemma}\label{sigmalemma}
The weak limit $\sigma$ is a solution of \eqref{eq4.9},
the block equation of the macroscopic model.
\end{lemma}

\begin{proof}
First take $\psi \in L^2(\Omega^3;
L^2(J;H^1_0(\Omega^{1,2}_m\times Q^3_m)))$ and define
\[
\bar{\psi}(x_1,x_2,x_3,z_3,t)
= \begin{cases}
\psi\big(x_1,x_2,x_3,\frac{z_3-c^{\varepsilon}(x_3)}{\varepsilon},t\big)
 & \text{if }  z_3\in \varepsilon Q^3_m+c^{\varepsilon}(x_3)\\
0 & \text{otherwise.}
\end{cases}
\]
Then, we have $\bar{\psi}\in L^2(\Omega^3;L^2(J;
H^1_0(\Omega^{\varepsilon}_m)))$.  We use $\bar{\psi}$
as the test function in the block equation of the weak form of
the microscopic model \eqref{eq5.2}.  So, for almost every
fixed $x_3\in \Omega^3$, we have
\[
(\phi^{\varepsilon}\sigma^{\varepsilon}_t,
\bar{\psi}(x_3))_{\Omega^{\varepsilon}_m\times J}
+(\widehat{\lambda}^{\varepsilon}\nabla_{x_1,x_2,z_3}
\sigma^{\varepsilon},\nabla_{x_1,x_2,z_3}
\bar{\psi}(x_3))_{\Omega^{\varepsilon}_m\times J}
= (f,\bar{\psi}(x_3))_{\Omega^{\varepsilon}_m\times J}.
\]
Since $\bar{\psi} =0$ in all blocks except the one that $x$ is in,
this simplifies to
\begin{align*}
&(\phi^{\varepsilon}\sigma^{\varepsilon}_t,
 \bar{\psi}(x_3))_{\Omega^{1,2}_m\times(\varepsilon Q^3_m
 + c^{\varepsilon}(x_3))\times J}\\
&+(\widehat{\lambda}^{\varepsilon}\nabla_{x_1,x_2,z_3}
 \sigma^{\varepsilon},\nabla_{x_1,x_2,z_3}
 \bar{\psi}(x_3))_{\Omega^{1,2}_m\times(\varepsilon Q^3_m
 + c^{\varepsilon}(x_3))\times J}\\
&= (f,\bar{\psi}(x_3))_{\Omega^{1,2}_m\times(\varepsilon
Q^3_m + c^{\varepsilon}(x_3))\times J}.
\end{align*}
Next, we let $y_3=\frac{z_3-c^{\varepsilon}(x_3)}{\varepsilon}$
and use the $ Q^3$-periodicity of $\phi$ and $\widehat{\lambda}$
 to obtain
\[
(\phi \widetilde{\sigma^{\varepsilon}_t},
\psi)_{\Omega^{1,2}_m\times Q^3_m\times J} +
(\widehat{\lambda}\nabla^{1/\varepsilon}_{x_1,x_2,y_3}
\widetilde{\sigma^{\varepsilon}},\nabla^{1/\varepsilon}
_{x_1,x_2,y_3}\psi)_{\Omega^{1,2}_m\times
Q^3_m\times J} = (\widetilde{f},\psi)_{\Omega^{1,2}_m\times
Q^3_m\times J},
\]
where $\nabla^{1/\varepsilon}_{x_1,x_2,y_3}
= (\frac{\partial}{\partial x_1},\frac{\partial}{\partial x_2},
\frac{1}{\varepsilon}\frac{\partial}{\partial y_3})$.
Simplifying and integrating in $x_3$, we get
\[
(\phi \widetilde{\sigma^{\varepsilon}_t}, \psi )_{\Omega_M\times
Q^3_m\times J} +(\lambda\nabla_{x_1,x_2,y_3}
\widetilde{\sigma^{\varepsilon}},\nabla_{x_1,x_2,y_3}
\psi)_{\Omega_M\times Q^3_m\times J} = (\widetilde{f},
\psi)_{\Omega_M\times Q^3_m\times J}.
\]
Finally, letting $\varepsilon \to 0$ gives us
\begin{equation}
(\phi \sigma_t, \psi )_{\Omega_M\times Q^3_m\times J}
+ (\lambda\nabla_{x_1,x_2,y_3}\sigma,\nabla_{x_1,x_2,y_3}
\psi)_{\Omega_M\times Q^3_m\times J}
= (f,\psi)_{\Omega_M\times Q^3_m\times J}.\label{eq10.20}
\end{equation}
This is a weak form of \eqref{eq4.9}.
\end{proof}

\begin{lemma}\label{rholemma}
The weak limit $\rho$ is a solution of \eqref{eq4.2}-\eqref{eq4.4},
\eqref{eq4.7}-\eqref{eq4.7a}, the fracture equations of
the macroscopic model.
\end{lemma}

\begin{proof}
Starting with \eqref{eq5.1}, we let $\varepsilon\to 0$ and consider
the convergence of each term.  The first and second terms are
straightforward:
\begin{gather*}
(\Phi^*\rho_t^{\varepsilon}, \varphi)_{\Omega^{\varepsilon}_f\times J}
= (\chi^{\varepsilon}_f\Phi^*\rho_t^{\varepsilon},
\varphi)_{\Omega\times J} \to (\Phi\rho_t, \varphi)_{\Omega\times J}, \\
(\Lambda\nabla\rho^{\varepsilon},
\nabla\varphi)_{\Omega^{\varepsilon}_f\times J}
=(\chi^{\varepsilon}_f\Lambda\nabla\rho^{\varepsilon},
\nabla\varphi)_{\Omega\times J} \to (\Lambda^{\#}\xi,
\nabla\varphi)_{\Omega\times J}.
\end{gather*}
Using the periodicity of $\phi$ and Lemmas \ref{lem7.3}-\ref{lem7.5},
the third term converges as follows:
\begin{align*}
(\phi^{\varepsilon}\sigma^{\varepsilon}_t,
\varphi)_{\Omega^{\varepsilon}_m\times J}
&= \frac{1}{|Q^3|}(\phi\widetilde{\sigma^{\varepsilon}_t},
\widetilde{\varphi})_{\Omega_M\times Q^3_m\times J} \\
& \to \frac{1}{|Q^3|}(\phi\sigma_t,\varphi)_{\Omega_M\times Q^3_m\times
J}
= -(f_m,\varphi)_{\Omega_M\times J}.
\end{align*}
Letting $\lambda_i^\varepsilon$, $i=1,2,3$, be the rows of
$\lambda^\varepsilon$, we expand the fourth term,
\begin{align*}
&(\widehat{\lambda}^{\varepsilon}\nabla\sigma^{\varepsilon},
\nabla\varphi)_{\Omega^{\varepsilon}_m\times J}\\
&=(\lambda_1^{\varepsilon}\cdot\nabla^\varepsilon
\sigma^{\varepsilon},\frac{\partial}{\partial
x_1}\varphi)_{\Omega^{\varepsilon}_m\times J} +
(\lambda_2^{\varepsilon}\cdot\nabla^\varepsilon\sigma^{\varepsilon},
\frac{\partial}{\partial x_2}\varphi)_{\Omega^{\varepsilon}_m\times J}
+ \varepsilon(\lambda_3^{\varepsilon}\cdot\nabla^\varepsilon
\sigma^{\varepsilon},\frac{\partial}{\partial x_3}\varphi)
_{\Omega^{\varepsilon}_m\times J}.
\end{align*}
The last term has an extra $\varepsilon$ and thus converges to $0$
by Lemma \ref{lem1}.  The convergence of other terms is shown
using Lemmas \ref{lem7.1}-\ref{lem7.5}:
\begin{align*}
(\lambda_i^{\varepsilon}\nabla^\varepsilon\sigma^{\varepsilon},
\frac{\partial}{\partial x_i}\varphi)_{\Omega^{\varepsilon}_m\times J}
&=\frac{1}{|Q^3|}(\widetilde{\lambda_i^{\varepsilon}}
 \widetilde{\nabla^\varepsilon\sigma^{\varepsilon}},
 \widetilde{\frac{\partial}{\partial x_i}\varphi})_{\Omega_M\times Q^3_m
 \times J}\\
&=\frac{1}{|Q^3|}(\lambda_{i}\cdot\nabla_{x_1,x_2,y_3}
\widetilde{\sigma^{\varepsilon}},
\widetilde{\frac{\partial}{\partial x_i}\varphi})_{\Omega_M\times Q^3_m
 \times J}\\
&\to \frac{1}{|Q^3|}(\lambda_i\cdot\nabla_{x_1,x_2,y_3}\sigma,
\frac{\partial}{\partial x_i}\varphi)_{\Omega_M\times Q^3_m\times J}\\
&= (f^i_M,\frac{\partial}{\partial
x_i}\varphi)_{\Omega_M\times J} \quad \text{for } i=1,2.
\end{align*}
Finally, we show the relationship between $\rho$ and $\xi$.
First, we have $\xi = \nabla\rho$ in $L^2(J;L^2(\Omega_F))$.
To see this, let $\varphi\in C^{\infty}_0(\Omega_F\times J)$
and extend it by $0$ to a function over all of $\Omega\times J$.
 Then for $i=1,2,3$, we have
\begin{align*}
0 &= (\chi_f^{\varepsilon}\frac{\partial}{\partial
x_i}\rho^{\varepsilon}, \varphi)_{\Omega\times J} -
(\chi_f^{\varepsilon}\frac{\partial}{\partial
x_i}\rho^{\varepsilon}, \varphi)_{\Omega_F\times J}\\
& \to (\frac{\Lambda^{\#}}{\Lambda}\xi_i, \varphi)_{\Omega\times J} -
(\frac{\Phi}{\Phi^*}\frac{\partial}{\partial x_i}\rho,
\varphi)_{\Omega_F\times J} \\
&= (\xi_i-\frac{\partial}{\partial x_i}\rho,
\varphi)_{\Omega_F\times J}.
\end{align*}
Similar calculations show that
$(\xi_1,\xi_2) = \nabla_{x_1,x_2}\rho
=(\frac{\partial}{\partial x_1}\rho,\frac{\partial}{\partial x_2}\rho)$
in $L^2(J;L^2(\Omega_M))$.  To relate $\xi_3$ and $\rho$ in
$\Omega_M\times J$, the following function will be useful.
Let $w\in C^0(\overline{ Q^3})\cap H^1( Q^3)$, with
$w\vert_{\overline{ Q^3_f}}\in C^1(\overline{ Q^3_f})$,
be a piecewise linear function that is $0$ on the boundary,
$\partial Q^3$, and constructed so that
$\frac{\partial}{\partial y_3}w=1$ in $\overline{ Q^3_f}$.
Then define a new function, $w^\varepsilon$, by
$w^{\varepsilon}(x_3) = \varepsilon w(\frac{x_3-c^{\varepsilon}
(x_3)}{\varepsilon})$ and note that
$w^\varepsilon\in H^1(\Omega^3)$. Also, note that for
$x_3\in\Omega^{3,\varepsilon}_f$, we have
$\frac{\partial}{\partial x_3}w^{\varepsilon}(x_3) = 1$.
Using Lemma \ref{lem7.3}, we see that
\[
\|w^{\varepsilon}\|_{L^2(\Omega\times J)}
= |Q^3|^{-1/2}\|\widetilde{w^{\varepsilon}}\|_{L^2(\Omega\times
Q^3\times J)}
= \varepsilon\Big(\frac{|\Omega|\,|J|}{|Q^3|}\Big)^{1/2}\|w\|_{L^2(
Q^3)} \to 0,
\]
and thus, $w^{\varepsilon}\to 0$ in $L^2(\Omega\times J)$.
Next, we let $\varphi \in C^{\infty}_0(\Omega_M\times J)$,
extend it by $0$ to all of $\Omega\times J$, and use
$w^{\varepsilon}\varphi$  as the test function in \eqref{eq5.1}:
\begin{align*}
&(\Phi^*\rho_t^{\varepsilon}, w^{\varepsilon}\varphi)_{\Omega^{\varepsilon}_f
\times J} + (\Lambda\nabla\rho^{\varepsilon}, \nabla(w^{\varepsilon}\varphi))
_{\Omega^{\varepsilon}_f\times J}\\
&+(\phi^{\varepsilon}\sigma^{\varepsilon}_t,
w^{\varepsilon}\varphi)_{\Omega^{\varepsilon}_m\times J} +
(\widehat{\lambda}^{\varepsilon}\nabla\sigma^{\varepsilon},
\nabla(w^{\varepsilon}\varphi))_{\Omega^{\varepsilon}_m\times
J} = (f,w^{\varepsilon}\varphi)_{\Omega\times J}.
\end{align*}
Letting $\varepsilon\to 0$, Lemma \ref{lem1} gives us
\begin{gather*}
(f,w^{\varepsilon}\varphi)_{\Omega\times J}\to 0,\quad
(\Phi^*\rho_t^{\varepsilon}, w^{\varepsilon}\varphi)_{
\Omega^{\varepsilon}_f\times J}\to 0,\\
(\phi^{\varepsilon}\sigma^{\varepsilon}_t,w^{\varepsilon}
\varphi)_{\Omega^{\varepsilon}_m\times J}\to 0,\quad
(\Lambda\nabla\rho^{\varepsilon},
w^{\varepsilon}\nabla\varphi)_{\Omega^{\varepsilon}_f\times
J}\to 0, \\
(\widehat{\lambda}^{\varepsilon}\nabla\sigma^{\varepsilon},
w^{\varepsilon}\nabla\varphi)_{\Omega^{\varepsilon}_m\times J}\to 0, \quad
(\widehat{\lambda}^{\varepsilon}\nabla\sigma^{\varepsilon},
\varphi\nabla w^{\varepsilon})_{\Omega^{\varepsilon}_m\times J}\to 0.
\end{gather*}
Also, because $\varphi = 0$ in $\Omega_F\times J$ and
$\frac{\partial}{\partial x_3}w^{\varepsilon} = 1$ in
$\Omega^{1,2}_m\times\Omega^{3,\varepsilon}_f\times J$, we have
\[
(\Lambda\nabla\rho^{\varepsilon}, \varphi\nabla
w^{\varepsilon})_{\Omega^{\varepsilon}_f\times J} =
(\Lambda\frac{\partial}{\partial x_3}\rho^{\varepsilon},
\varphi\frac{\partial}{\partial
x_3}w^{\varepsilon})_{\Omega^{\varepsilon}_f\times J}=
(\Lambda\frac{\partial}{\partial x_3}\rho^{\varepsilon},
\varphi)_{\Omega^{1,2}_m\times\Omega^{3,\varepsilon}_f\times J}.
\]
Combining the above gives us
$(\Lambda\frac{\partial}{\partial x_3}\rho^{\varepsilon},
\varphi)_{\Omega^{1,2}_m\times\Omega^{3,\varepsilon}_f\times J} \to 0$.
But,
\[
(\Lambda\frac{\partial}{\partial x_3}\rho^{\varepsilon},
\varphi)_{\Omega^{1,2}_m\times\Omega^{3,\varepsilon}_f\times J} =
(\chi^{\varepsilon}_f\Lambda\frac{\partial}{\partial
x_3}\rho^{\varepsilon}, \varphi)_{\Omega\times J} \to
(\Lambda^\#\xi_3, \varphi)_{\Omega\times J}=(\Lambda_M\xi_3,
\varphi)_{\Omega_M\times J}.
\]
Thus, $(\xi_3, \varphi)_{\Omega_M\times J} = 0$ for every
$\varphi\in C^{\infty}_0(\Omega_M\times J)$ and $\xi_3=0$
in $\Omega_M\times J$.  Putting all the convergence results together,
we see that
\begin{equation}
\begin{aligned}
&(\Phi\rho_t, \varphi)_{\Omega\times J} + (\Lambda_F\nabla\rho, \nabla\varphi)
_{\Omega_F\times J} + (\Lambda_M\nabla_{x_1,x_2}\rho, \nabla_{x_1,x_2}
\varphi)_{\Omega_M\times J}\\
&+ (f^1_M,\frac{\partial}{\partial x_1}\varphi)_{\Omega_M\times J}
  + (f^2_M,\frac{\partial}{\partial x_2}\varphi)_{\Omega_M\times J}\\
&= (f,\varphi)_{\Omega\times J} + (f_m,\varphi)_{\Omega_M\times J},
\quad \forall\ \varphi\in L^2(J;H^1(\Omega)),
\end{aligned}\label{eq11.44}
\end{equation}
which is a weak form of \eqref{eq4.2}-\eqref{eq4.4},
\eqref{eq4.7}-\eqref{eq4.7a}.
\end{proof}

Next we show that $\rho$ and $\sigma$ satisfy the appropriate
initial and boundary conditions.  Throughout we will
use the following well-known results.

\begin{proposition}\label{lem13.2}
If $X$ and $Y$ are Banach spaces, $T:X\to Y$ a bounded linear
operator, and $\psi_{\alpha}\rightharpoonup \psi$ weakly in $X$,
then $T(\psi_{\alpha})\rightharpoonup T(\psi)$ weakly in $Y$.
\end{proposition}

\begin{proposition}\label{lem13.21}
Let $U$ be a bounded domain (not sitting on both sides of
$\partial U$) and $U_1, U_2$ be disjoint open subsets of $U$
with $\partial U_1\cap\partial U_2=\Sigma$ ($\Sigma\neq\emptyset$
 and Lipschitz continuous).  Furthermore, let $\Gamma$ be any
subset of $\Sigma$, including $\Sigma$ itself, and $u\in H^1(U)$.
If $T_{U_1}:H^1(U_1)\to L^2(\Gamma)$, and
$T_{U_2}:H^1(U_2)\to L^2(\Gamma)$ are trace operators, then
$T_{U_1}(u\vert_{U_1})=T_{U_2}(u\vert_{U_2})$ in $L^2(\Gamma)$.
\end{proposition}

We first show the initial conditions.

\begin{lemma}\label{initcond}
The weak limits $\rho$ and $\sigma$ satisfy the initial conditions,
\eqref{eq4.6} and \eqref{eq4.11}, of the macroscopic model:
$\rho(0)=\rho_{\rm init}$ in $L^2(\Omega)$ and
$\sigma(0)=\rho_{\rm init}$ in $L^2(\Omega^3;L^2(\Omega^{1,2}_m\times
Q^3_m))$.
\end{lemma}

\begin{proof}
There exist appropriate time-zero trace operators,
$T^f_0:H^1(J;L^2(\Omega))\to L^2(\Omega)$ and
$T^m_0:L^2(\Omega^3;H^1(J; L^2(\Omega^{1,2}_m\times Q^3_m)))
\to L^2(\Omega^3;L^2(\Omega^{1,2}_m\times Q^3_m))$, that are
linear and bounded.  Thus, \eqref{eq9.1} and
Proposition \ref{lem13.2} give us
\[
T^f_0(\chi^{\varepsilon}_f \Phi^*\rho^{\varepsilon})
\rightharpoonup T^f_0(\Phi \rho) = \Phi T^f_0\rho \quad
\text{weakly in}\ L^2(\Omega).
\]
Also, recalling the initial condition \eqref{eq5.4}, we have
\[
T^f_0(\chi^{\varepsilon}_f \Phi^*\rho^{\varepsilon})
=\chi^{\varepsilon}_f
\Phi^*\rho^{\varepsilon}(0)=\chi^{\varepsilon}_f
\Phi^*\rho_{\rm init}\rightharpoonup \Phi \rho_{\rm init} \quad
\text{weakly in } L^2(\Omega).
\]
Therefore, $T^f_0\rho=\rho_{\rm init}$ in $L^2(\Omega)$.
 Similarly, by \eqref{eq9.7} and Proposition \ref{lem13.2} we have
\[
T^m_0\widetilde{\sigma}^{\varepsilon} \rightharpoonup
T^m_0\sigma\quad\text{weakly in }
L^2(\Omega^3;L^2(\Omega^{1,2}_m\times Q^3_m)).
\]
But the initial condition \eqref{eq5.4} and Lemma \ref{lem7.5}
give us
\[
T^m_0\widetilde{\sigma}^{\varepsilon}
=\widetilde{\sigma}^{\varepsilon}(0)=\widetilde{\rho_{\rm init}}\to
\rho_{\rm init}\quad\text{strongly in }
L^2(\Omega^3;L^2(\Omega^{1,2}_m\times Q^3_m)).
\]
Hence, $T^m_0\sigma=\rho_{\rm init}$ in
$L^2(\Omega^3;L^2(\Omega^{1,2}_m\times Q^3_m))$.
\end{proof}

It remains to show \eqref{eq4.10}.
To do so, we must consider the boundary of the block in
pieces because the domain of $\widetilde{\rho^\varepsilon}$
depends on $x_1$ and $x_2$.  We start with the lateral boundaries
of the block, $\partial\Omega^{1,2}_m\times\overline{Q^3_m}$.

\begin{lemma}\label{lemlateral}
The weak limits $\rho$ and $\sigma$ satisfy $\rho =\sigma$
in $L^2(\Omega^3;L^2(\partial\Omega^{1,2}_m\times\overline{Q^3_m}
\times J))$.
\end{lemma}

\begin{proof}
Let $T^{a,\varepsilon}_{f_1}$ and  $T^{a,\varepsilon}_{m_1}$ be
standard traces for the fracture and matrix domains respectively:
$T^{a,\varepsilon}_{f_1}:H^1(\Omega^{\varepsilon}_f\times J)
\to L^2(\partial\Omega^{1,2}_m\times\overline{\Omega_m^{3,
\varepsilon}}\times J)$ and $T^{a,\varepsilon}_{m_1}
:H^1(\Omega^{\varepsilon}_m\times J)\to L^2(\partial
\Omega^{1,2}_m\times\overline{\Omega_m^{3,\varepsilon}}\times J)$.

Additionally, let $T^a_f:L^2(\Omega^3;H^1(\Omega^{1,2}_f\times
Q^3\times J))\to L^2(\Omega^3;L^2(\partial\Omega^{1,2}_m
\times\overline{Q^3_m}\times J))$ and $T^a_m:L^2(\Omega^3;
H^1(\Omega^{1,2}_m\times Q^3_m\times J))\to L^2(\Omega^3;L^2
(\partial\Omega^{1,2}_m\times\overline{Q^3_m}\times J))$
be standard trace operators. Then the dilation and trace
operators commute in the following sense:
\begin{equation}
T^a_f \widetilde{\rho}^{\varepsilon}
= \widetilde{T^{a,\varepsilon}_{f_1} \rho^{\varepsilon}},\quad
 T^a_m \widetilde{\sigma}^{\varepsilon}
= \widetilde{T^{a,\varepsilon}_{m_1} \sigma^{\varepsilon}}\quad
\text{in } L^2(\Omega^3;L^2(\partial\Omega^{1,2}_m\times
\overline{Q^3_m}\times J))\label{result1} .
\end{equation}
This result is straightforward for all functions in
$C^0(\overline{\Omega_r^{\varepsilon}\times J})
\cap H^1(\Omega_r^{\varepsilon}\times J)$, $r=f$ or $m$.
The full result then follows from a density argument since
$C^0(\overline{\Omega_r^{\varepsilon}\times J})$ is dense
in $H^1(\Omega_r^{\varepsilon}\times J)$.  Next, we use
Proposition \ref{lem13.21} with $U=\Omega\times J$,
$U_1=\Omega^{\varepsilon}_f\times J$,
$U_2=\Omega^{\varepsilon}_m\times J$, and
$\Gamma=\partial\Omega^{1,2}_m\times\overline{\Omega^{3,\varepsilon}_m}
\times J$, to get
\[
T^{a,\varepsilon}_{f_1} \rho^{\varepsilon}=T^{a,\varepsilon}_{m_1}
\sigma^{\varepsilon}\quad\text{in }
L^2(\partial\Omega^{1,2}_m\times\overline{\Omega^{3,\varepsilon}_m}\times
J).
\]
Then dilating both sides, we have
\begin{equation}
\widetilde{T^{a,\varepsilon}_{f_1} \rho^{\varepsilon}}
=\widetilde{T^{a,\varepsilon}_{m_1} \sigma^{\varepsilon}}\quad
\text{in }  L^2(\Omega^3;L^2(\partial\Omega^{1,2}_m\times
\overline{Q^3_m}\times J))\label{result2}.
\end{equation}
Together, \eqref{result1} and \eqref{result2} imply that
$T^a_f \widetilde{\rho}^{\varepsilon}
=T^a_m \widetilde{\sigma}^{\varepsilon}$ in
$L^2(\Omega^3;L^2(\partial\Omega^{1,2}_m\times\overline{Q^3_m}\times J))$.
 Since the trace operators are all linear and bounded,
Proposition \ref{lem13.2} gives us
\begin{gather*}
T^a_f \widetilde{\rho}^{\varepsilon}\rightharpoonup
T^a_f\tau_3 = T^a_f\rho \quad\text{weakly in }
 L^2(\Omega^3;L^2(\partial\Omega^{1,2}_m\times\overline{Q^3_m}\times J)),\\
T^a_m \widetilde{\sigma}^{\varepsilon}\rightharpoonup
T^a_m\sigma \quad\text{weakly in }
L^2(\Omega^3;L^2(\partial\Omega^{1,2}_m\times\overline{Q^3_m}\times
J)),
\end{gather*}
and the main result follows.
\end{proof}

Next we derive a similar boundary condition for the top and
bottom of the blocks, $\overline{\Omega^{1,2}_m}\times\partial Q^3_m$.
The proof is nearly identical to the one for the lateral boundaries.

\begin{lemma}\label{lemtopbot}
The weak limits $\rho$ and $\sigma$ satisfy $\rho =\sigma$
in $L^2(\Omega^3;L^2(\overline{\Omega^{1,2}_m}\times\partial
Q^3_m\times J))$.
\end{lemma}

\begin{proof}
In similar fashion, we let $T^{b,\varepsilon}_{f_1}:
H^1(\Omega^{\varepsilon}_f\times J)\to L^2(\overline{\Omega^{1,2}_m}
\times\partial\Omega_m^{3,\varepsilon}\times J)$,
$T^{b,\varepsilon}_{m_1}:H^1(\Omega^{\varepsilon}_m\times J)
\to L^2(\overline{\Omega^{1,2}_m}\times\partial
\Omega_m^{3,\varepsilon}\times J)$,
$T^b_f:L^2(\Omega^3;H^1(\Omega^{1,2}_m\times Q^3_f\times J))
\to L^2(\Omega^3;L^2(\overline{\Omega^{1,2}_m}\times\partial
Q^3_m\times J))$ and  $T^b_m:L^2(\Omega^3;H^1(\Omega^{1,2}_m\times
Q^3_m\times J))\to L^2(\Omega^3;L^2(\overline{\Omega^{1,2}_m}
\times\partial Q^3_m\times J))$ be standard trace operators.
As before, we need the dilation and trace operators to commute:
\begin{equation}
T^b_f \widetilde{\rho}^{\varepsilon}
= \widetilde{T^{b,\varepsilon}_{f_1} \rho^{\varepsilon}},\quad
 T^b_m \widetilde{\sigma}^{\varepsilon}
= \widetilde{T^{b,\varepsilon}_{m_1} \sigma^{\varepsilon}}
\quad\text{in } L^2(\Omega^3;L^2(\overline{\Omega^{1,2}_m}\times\partial
Q^3_m\times J))\label{result3}.
\end{equation}
Again, the result holds for all functions in
$C^0(\overline{\Omega_r^{\varepsilon}\times J})\cap
H^1(\Omega_r^{\varepsilon}\times J)$, $r=f$ or $m$,
and thus for all functions in $H^1(\Omega_r^{\varepsilon}\times J)$
via a density argument.  Next, we use Proposition \ref{lem13.21}
with $U=\Omega\times J$, $U_1=\Omega^{\varepsilon}_f\times J$,
$U_2=\Omega^{\varepsilon}_m\times J$, and
$\Gamma=\overline{\Omega^{1,2}_m}\times\partial
\Omega^{3,\varepsilon}_m\times J$ to get
\[
T^{b,\varepsilon}_{f_1} \rho^{\varepsilon}=T^{b,\varepsilon}_{m_1}
\sigma^{\varepsilon},\quad
L^2(\overline{\Omega^{1,2}_m}\times\partial\Omega^{3,\varepsilon}_m\times
J).
\]
Thus, dilating both sides gives
\begin{equation}
\widetilde{T^{b,\varepsilon}_{f_1} \rho^{\varepsilon}}
=\widetilde{T^{b,\varepsilon}_{m_1} \sigma^{\varepsilon}}
\quad\text{in } L^2(\Omega^3;L^2(\overline{\Omega^{1,2}_m}
\times\partial Q^3_m\times J))\label{result4}.
\end{equation}
Finally, \eqref{result3} and \eqref{result4} together imply
that $T^b_f \widetilde{\rho}^{\varepsilon}=T^b_m
\widetilde{\sigma}^{\varepsilon}$ in $L^2(\Omega^3;
L^2(\overline{\Omega^{1,2}_m}\times\partial Q^3_m\times J))$.
Then the main result follows from Proposition \ref{lem13.2} since
\begin{gather*}
T^b_f \widetilde{\rho}^{\varepsilon}\rightharpoonup T^b_f\tau
= T^b_f\rho \quad\text{weakly in}\ L^2(\Omega^3;
L^2(\overline{\Omega^{1,2}_m}\times\partial Q^3_m\times J)),\\
T^b_m \widetilde{\sigma}^{\varepsilon}\rightharpoonup
T^b_m\sigma \quad\text{weakly in }
L^2(\Omega^3;L^2(\overline{\Omega^{1,2}_m}\times\partial
Q^3_m\times J)).
\end{gather*}
\end{proof}

Thus, from Lemmas \ref{lemlateral} and \ref{lemtopbot},
the limit solution satisfies the full boundary condition,
\eqref{eq4.10}.  We now show that our limit solution is unique.

\begin{lemma}\label{thm13.1}
The weak solutions $\rho$ and $\sigma$ are unique.
\end{lemma}

\begin{proof}
Let $\rho$ and $\sigma$ be the difference of two solutions to
macroscopic model.  Then, $\rho$ and $\sigma$ satisfy
\eqref{eq10.20} and \eqref{eq11.44} with the source terms involving
$f$ equal to zero.  Additionally, $\rho$ and $\sigma$ satisfy
\eqref{eq4.10}, $\rho(0)=0$, and $\sigma(0)=0$. Next, using
$\rho$ and $\sigma-\rho$ as test functions, we see that $\rho$
and $\sigma$ satisfy the following:
\begin{equation}
\begin{aligned}
&(\Phi\rho_t, \rho)_{\Omega\times J}
+ (\Lambda_F\nabla\rho, \nabla\rho)_{\Omega_F\times J}
+ (\Lambda_M\nabla_{x_1,x_2}\rho, \nabla_{x_1,x_2}
 \rho)_{\Omega_M\times J}\\
&+ (f^1_M,\frac{\partial}{\partial x_1}\rho)_{\Omega_M\times J}
 + (f^2_M,\frac{\partial}{\partial x_2}\rho)_{\Omega_M\times J}
 = (f_m,\rho)_{\Omega_M\times J}\label{uniq1}
\end{aligned}
\end{equation}
and
\begin{equation}
\begin{aligned}
&(\phi \sigma_t, \sigma)_{\Omega_M\times Q^3_m\times J}
-(\phi \sigma_t, \rho )_{\Omega_M\times Q^3_m\times J}
+ (\lambda\nabla_{x_1,x_2,y_3}\sigma,\nabla_{x_1,x_2,y_3}
 \sigma)_{\Omega_M\times Q^3_m\times J}\\
&-(\lambda\nabla_{x_1,x_2,y_3}\sigma,\nabla_{x_1,x_2,y_3}
 \rho)_{\Omega_M\times Q^3_m\times J}=0.
\end{aligned} \label{uniq2}
\end{equation}
Then multiplying \eqref{uniq1} by $|Q^3|$, adding it to
\eqref{uniq2}, and recalling what $f_m$, $f^1_M$, and $f^2_M$ are,
gives us
\begin{align*}
&|Q^3|(\Phi\rho_t,\rho)_{\Omega\times J} + (\phi\sigma_t,\sigma)_{\Omega_M
\times Q^3_m\times J} + |Q^3|(\Lambda_F\nabla\rho,\nabla\rho)_{\Omega_F\times J} \\
&+ |Q^3|(\Lambda_M\nabla_{x_1,x_2}\rho,\nabla_{x_1,x_2}
\rho)_{\Omega_M\times J} +
(\lambda\nabla_{x_1,x_2,y_3}\sigma,\nabla_{x_1,x_2,y_3}
\sigma)_{\Omega_M\times Q^3_m\times J} = 0 .
\end{align*}
Hence, standard energy estimates of this equation give us
$\|\rho\|_{L^2(\Omega\times J)} = 0$ and
$\|\sigma\|^2_{L^2(\Omega_M\times Q^3_m\times J)} = 0$, and the
result follows.
\end{proof}

Since the solution is unique, the whole sequence of solutions
to the microscopic model converges to $\rho$ and $\sigma$.
Finally, we address the well-posedness of the macroscopic model.

\begin{theorem}\label{thm9.0}
The weak limits $\rho$ and $\sigma$ depend continuously on the data.
That is, $\|\rho\|_{H^1(J;L^2(\Omega))}
+ \|\sigma\|_{L^2(\Omega^2;H^1(\Omega^1_m\times Q^3_m\times J))}
\leq C \left(\|f\|_{L^2(J;L^2(\Omega))}
+ \|\rho_{\rm init}\|_{L^2(\Omega)}\right)$.
\end{theorem}

\begin{proof}
It follows from the energy estimates used to prove Theorem \ref{thm1}
that $\rho^\varepsilon$ and $\sigma^\varepsilon$ depend continuously
on the data, $f$ and $\rho_{\rm init}$.  An appropriate bound on
the norm of $\widetilde{\sigma^\varepsilon}$ can be obtained by
performing calculations similar to those found in the proof of
Lemma \ref{lem7.3}. The result then follows by passing to the
weak limits.
\end{proof}

\section{Extensions}\label{sec6}

The geometry considered in this work can be extended to include
larger fractures in the vertical direction in at least two ways.
 First of all, the domain of consideration would consist of
horizontal translations of cells, where each cell is a layered
porous medium of the type described in this work (See Fig. \ref{fig15}).

\begin{figure}[ht]
  \begin{center}
  \includegraphics[scale=.4]{fig6.jpg}
  \end{center}
\caption{A vertical cross-section of a domain with large vertical fractures 
in between horizontal translations of cells with the same structure as 
in Fig. \ref{fig1}.}\label{fig15}
\end{figure}

One method of homogenizing this medium would be to use the following
recursive homogenization approach:  first homogenize each cell in
the vertical direction using the method in this work; then homogenize
the vertical fractures using a homogenization approach in the
horizontal direction on a larger length scale.

A second approach would be to use reiterative homogenization.
That is, we would set up two different length scales, both dependent
on a parameter $\varepsilon$, that would eventually tend to zero in
the homogenization process.  The fine scale would be the horizontal
fractures in each cell and the coarse scale would be the larger
vertical fractures between the cells.  As $\varepsilon\to 0$, the
fine scale fractures would tend to zero faster than the coarse
scale fractures.

\subsection*{Acknowledgements}
  The authors would like to express their sincere gratitude
to Steve Wright for his enthusiastic conversations on homogenization
techniques.  They would also like to thank the reviewers of the
manuscript for their thoughtful suggestions that led to a better
presentation.

\begin{thebibliography}{00}

\bibitem{ArDH90} Arbogast, T.; Douglas Jr., J.; Hornung, U.;
 Derivation of the double porosity model of single phase flow via
homogenization theory,
\emph{SIAM J. Math. Anal.}, {\bf 21} (1990), 823--836.

\bibitem{ArDoHu91} Arbogast, T.; Douglas Jr., J.;  Hornung, U.;
Modeling of naturally fractured reservoirs by formal homogenization
techniques, \emph{Frontiers in Pure and Applied Mathematics},
R. Dautray (ed.), Elsevier, Amsterdam, 1991,  1--19.

\bibitem{Ar89} Arbogast, T.;
Analysis of the simulation of single phase flow through a naturally
fractured reservoir, \emph{SIAM J. Numer. Anal.},
{\bf 26} (1989), 12--29.

\bibitem{BZK60} Barenblatt, G. I.; Zheltov, P.;  Kochina, I. N.;
 Basic concepts in the theory of seepage of homogeneous liquids
in fissured rocks,
\emph{Prikl. Mat. Mekh.}, {\bf 24} (1960), 852--864
(in Russian); \emph{J. Appl. Math. Mech.}, {\bf 24} (1960), 1286--1303.


\bibitem{Be72} Bear, J.;
 \emph{Dynamic of Fluids in Porous Media,} Dover Publications,
New York, 1972.

\bibitem{ClSh} Clark, G.; Showalter, R. E.;
 Fluid flow in layered media,
\emph{Quart. Appl. Math.}, {\bf 52} (1994), no. 4, 777--795.

\bibitem{da90} Douglas Jr. J.; Arbogast, T.;
Dual porosity models for flow in naturally fractured reservoirs,
in \emph{Dynamics of Fluids in Hierarchical Porous Formations},
J.\, H. Cushman, ed., Academic Press, London, 1990, 177--221.

\bibitem{DoAbPLHeNu} Douglas Jr., J.; Arbogast, T.; Paes-Leme, P. J.;
 Hensley, J.; Nunes, N.;
 Immiscible displacement in vertically fractured reservoirs,
\emph{Transport in Porous Media}, {\bf 12} (1993), 73--106.

\bibitem{DoPeSh} Douglas Jr., J.; Peszy\'{n}ska, M.; Showalter, R. E.;
Single phase flow in partially fissured media,
\emph{Transport in Porous Media}, {\bf 28} (1997), 285--306.

\bibitem{DPLAS87} Douglas Jr., J.; Paes Leme, P. J.;
 Arbogast T.; Schmitt, T.;
 Simulation of flow in naturally fractured reservoirs,
\emph{Proc. Ninth SPE Symposium on Reservoir Simulation},
Dallas, 1987, 271--279.

\bibitem{DKPLS98} Douglas Jr., J.; Kischinhevsky, M.; Paes Leme, P. J.;
Spagnuolo, A.;
 A multiple-porosity model for a single-phase flow through
naturally-fractured porous media, \emph{Comp. Appl. Math.},
{\bf 17} (1998), 19--48.

\bibitem{HoBook} Hornung, U.;
 \emph{Homogenization and Porous Media}, Springer, New York, 1997.

\bibitem{K69} Kazemi, H.;
 Pressure transient analysis of naturally fractured reservoirs
with uniform fracture distribution,
\emph{J. Soc. Petroleum Engrs.}, {\bf 9} (1969), 451--462.

\bibitem{Pirson} Pirson, S. J.;
 Performance of fractured oil reservoirs,
\emph{Bull. Amer. Assoc. Petroleum Geologists},
\textbf{37} (1953), pp. 232-244.

\bibitem{SaPa} Sanchez-Palencia, E.;
 Non-homogeneous Media and Vibration Theory,
\emph{Lecture Notes in Physics}, 127, Springer-Verlag, New York, 1980

\bibitem{ShSpWr} Shi, P.; Spagnuolo, A.;  Wright, S.;
Reiterated homogenization and the double-porosity model,
\emph{Transport in Porous Media}, {\bf 59} (2005), 73--95.


\bibitem{SpW1} Spagnuolo, A. M.;  Wright, S.;
 Analysis of a multiple-porosity model for single-phase
flow through naturally-fractured porous media, \emph{J. Appl. Math},
{\bf 7} (2003), 327--364.

\bibitem{SpW2} Spagnuolo, A. M.; Wright, S.;
 Derivation of a multiple porosity model of single-phase flow
through a fractured porous medium via recursive homogenization,
 \emph{Asymptotic Anal.}, {\bf 39} (2004), 91--112.

\bibitem{WR63} Warren, J. E. and Root, P. J.;
 The behavior of naturally fractured reservoirs,
\emph{J. Soc. Petroleum Engrs.}, {\bf 3} (1963), 245--255.

\end{thebibliography}

\end{document}