\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: \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} 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 \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} 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 $$(\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}$$ 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{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} 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: $$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} .$$ 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 $$\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}.$$ 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: $$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}.$$ 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 $$\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}.$$ 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{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} and \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} 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}