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

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

\begin{document}
\title[\hfilneg EJDE-2014/189\hfil Geological fissured systems]
{Homogenization of geological fissured systems with curved non-periodic cracks}

\author[F. A. Morales\hfil EJDE-2014/189\hfilneg]
{Fernando A. Morales}  % in alphabetical order

\address{Fernando A. Morales \newline
Escuela de Matem\'aticas, Universidad Nacional de Colombia,
Sede Medell\'in, Colombia}
\email{famoralesj@unal.edu.co}


\thanks{Submitted September 7, 2013. Published September 11, 2014.}
\subjclass[2000]{35F15, 80M40, 76S99, 35B25}
\keywords{Fissured media; tangential flow; interface geometry;
\hfill\break\indent  coupled Darcy flow system; upscaling; mixed formulation}

\begin{abstract}
 We analyze the steady fluid flow in a porous medium containing a network
 of thin fissures of width $\mathcal{O}(\epsilon)$, generated by the rigid
 translation of continuous piecewise $C^{1}$ functions in a fixed direction.
 The phenomenon is modeled in mixed variational formulation, using the
 stationary Darcy's law and coefficients of low resistance $\mathcal{O}(\epsilon)$
 on the network. The singularities are removed by asymptotic analysis as
 $\epsilon \to 0$ which yields an analogous system hosting only tangential
 flow in the fissures. Finally the fissures are collapsed into two dimensional
 manifolds.
\end{abstract}

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

\section{Introduction}

Groundwater and oil reservoirs are frequently fissured or layered; i.e.,
the bedrock contains fissures of characteristic dimensions considerably higher
than those of the average pore size of the rock. The modeling of saturated
flow through geological structures such as these, gives rise to singular
problems of partial differential equations \cite{Show97}. On one hand,
the singularities are caused by the drastic change of permeability from the
rock matrix to the fissures. On the other hand, a geometric singularity
is introduced due to the thinness of the fractures. The presence of singularities
in the model has non-desirable effects in their numerical implementation;
some of these are ill-conditioned matrices, high computational costs,
numerical instability, etc. This subject is a very active research field,
see \cite{ArbBrunson2007, Gatica2009, Yotov, Loredana1, JaffRob05}
for numerical analysis aspects, see \cite{ZhaoQin, Levy83} for
modeling discussion and see
\cite{Allaire2009, ArbLehr2006, Gunzburger2009, Mikelic89, Morales1, MoralesShow1}
for rigorous mathematical treatment of the phenomenon. Homogenization and
asymptotic analysis techniques are a common approach for the analytical point of view.
However, the remarkable achievements in the field require very restrictive
hypotheses for the description of the geometry such as regular geometric shapes
or periodically arrayed structures \cite{Hornung, SP80}.
In general the variational methods for partial differential equations can give
successful formulations for a wide class of geometric domains.
The limited treatment of the geometry is due to the notorious technical difficulties,
that more general shapes introduce in the asymptotic analysis of the problem.

In the present work, the geometric possibilities of the medium are extended
to fissures that are not necessarily flat. We use the
\emph{mixed mixed formulation} and the scaling for the flow resistance
coefficients presented in \cite{MoralesShow2}, then a careful choice of
directions or ``stream lines'', consistent with the natural geometry of the
problem permits a successful asymptotic analysis of the model.
This leads to a system coupled through multiple two dimensional manifolds
representing the fissures in the upscaled model. Additionally, the formulation
allows remarkable generality in the fluid exchange balance conditions between
the rock matrix and the channels and substantial efficiency for handling the
system of equations as well as the information (coefficients, matrices, etc)
describing the geometry of the fractures. This is mostly due to the fact that
the formulation does not demand coupling constraints on the underlying spaces
of functions. The main goal of the paper is to emphasize the geometry,
consequently the study is limited to the steady case. We describe flow
with Darcy's law
\begin{subequations}\label{Pblm FORMAL fissures system strong problem}
\begin{equation}\label{Eq FORMAL Darcy}
a (\cdot) \mathbf{u} + \boldsymbol{\nabla} p + \mathbf{g} = 0 ,
\end{equation}
together with the conservation law
\begin{equation}\label{Eq FORMAL fluid mass conservation}
\boldsymbol{\nabla}\cdot \mathbf{u} = F.
\end{equation}
Drained and null-flux boundary conditions will be specified on different
 parts of the boundary of the domain to set a boundary value problem. 
The balance conditions of normal stress and normal flux across the 
interface $\Gamma$, separating the regions $\Omega_{1}$ and $\Omega_{2}$ 
are given by
\begin{gather}\label{Eq FORMAL interface balance 1 gamma}
 p^1 -  p^2 = \alpha   \mathbf{u}^1 \cdot \boldsymbol{\widehat{n}}\quad \text{and}\\
\label{Eq FORMAL interface balance 2 gamma}
\mathbf{u}^1 \cdot \boldsymbol{\widehat{n}} 
- \mathbf{u}^2\cdot\boldsymbol{\widehat{n}}  =
f_{\Gamma} \quad \text{on }  \Gamma .
\end{gather}
\end{subequations}
The superscripts denote restrictions to $\Omega_{1}$ and $\Omega_{2}$.
The coefficient $a(\cdot)$ is the flow resistance, i.e. the fluid viscosity
times the inverse of the permeability of the medium, to be scaled consistently
with the fast and slow flow regions of the medium. Finally, the coefficient
$\alpha$ indicates the fluid entry resistance of the rock matrix.

In the following section we define the geometric setting, formulate the problem
in mixed mixed variational formulation and prove its well-posedness.
In section three the problem is referred to a common geometric setting in
order make possible the asymptotic analysis, the existence of a-priori estimates
and the structure of the limiting solution are also shown. Section four studies
the formulation and well-posedness of the limiting problem as well as its strong
form (particularly important for boundary and interface conditions); it also
proves the strong convergence of the solutions. Section five sets the limiting
problem as a coupled system with two dimensional interfaces and section six
discusses the possibilities and limitations of the technique as well as related
future work.

\section{Formulation and geometric setting}

Vectors are denoted by boldface letters as are vector-valued functions
and corresponding function spaces. We use $\widetilde{\mathbf{x}}$ to indicate a vector
in $\mathbb{R}^2$; if $\mathbf{x}\in \mathbb{R}^3$ then the $\mathbb{R}^2\times\{0\}$
projection is identified with $\widetilde{\mathbf{x}}:=(x_{1}, x_{2})$ so that
$\mathbf{x} = (\widetilde{\mathbf{x}},  x_3)$. The symbol $\boldsymbol{\widetilde{\nabla}}$ represents the gradient in
the first first two directions: $\boldsymbol{\widehat{\imath}}$, $\boldsymbol{\widehat{\jmath}}$. Given a function
$f: \mathbb{R}^3\to \mathbb{R}$ then  $\int_{\mathcal {M} } f dS$ is
the notation for its surface integral on the two dimensional manifold
$\mathcal{M}\subseteq \mathbb{R}^3$. $\int_{A} f\,d\mathbf{x}$ stands for the volume
integral in the set $A\subseteq \mathbb{R}^3$; whenever the context is clear
we simply write $\int_{A} f$.

The symbol $\boldsymbol{\widehat{\nu} }$ denotes the outwards normal vector
on the boundary of a given domain $\mathcal{O}$ and $\boldsymbol{\widehat{n}}$ denotes the upwards
vector normal to a given surface i.e. $\boldsymbol{\widehat{n}} \cdot \boldsymbol{\widehat{k}} \geq 0$.
For any $A\subseteq\mathbb{R}^{3}$ and $t\in \mathbb{R}$ we define its $t$-vertical shift by
\begin{equation}\label{Def general vertical shift}
A+t:=\big\{\mathbf{x}+ t \boldsymbol{\widehat{k}}: \mathbf{x}\in A\big\}
\end{equation}

\subsection{General geometric setting}

The present work will be limited to the study of fractured media where each
fissure can be described in a specific way.

\begin{definition}\label{Def eligible surface and generated fissure} \rm
Let $\mathcal{O}\subseteq \mathbb{R}^{2}$ be a bounded open simply connected
set and $\zeta\in C(\overline{\mathcal{O}})$ be a piecewise $C^{1}$ function.
Define the surface
\begin{equation}\label{Def C2 surface}
\Upsilon:= \{[\widetilde{\mathbf{x}}, \zeta(\widetilde{\mathbf{x}})]
:\widetilde{\mathbf{x}}\in \mathcal{O}\} .
\end{equation}
We assume that $\Upsilon$ is non-vertical with 
$\operatorname{ess\,inf}\{ \boldsymbol{\widehat{n}}(s)\cdot\boldsymbol{\widehat{k}}:
 s\in \Upsilon\}>0$. Given $h>0$, define the fissure of height $h$ generated 
by a rigid vertical translation of $\Upsilon$ as the domain
\begin{equation}\label{Def vertical fissure}
\Theta(h, \Upsilon):= \{(\widetilde{\mathbf{x}}, y):
\zeta(\widetilde{\mathbf{x}})< y < \zeta(\widetilde{\mathbf{x}})+h\} .
\end{equation}
\end{definition}


\begin{remark}\label{Rem comments on the height width relation} \rm
Notice that in the definition of $\Theta(h, \Upsilon)$ we mention $h$ as
the height and not as the width of the crack.
Figure \ref{Fig Fissure with high gradients} shows that, depending on
the gradient of the surface, the height $h$ can become significantly
different from the actual width.
\end{remark}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig1}
\end{center} % UnidirectionalSystem3.pdf
\caption{Unidirectional translation generated fissures}\label{Fig Fractured Medimu}
\end{figure}

Figure \ref{Fig Fractured Medimu} depicts a region
$\Omega\subseteq\mathbb{R}^{3}$ containing a network of fissures generated
by vertical rigid translation of continuous piecewise $C^{1}$ surfaces.
For the sake of clarity we restrict the analysis to the case where only one
fissure is embedded in the domain. From here, the generalization to a system
with multiple fissures is a rather simple exercise. Next we define this domain

\begin{definition}\label{Def unidirectional translation medium} \rm
A medium $\Omega$ with a vertical translation generated fissure is composed
by three connected components denoted $\Theta_{1}$, $\Theta_{2}$ and
$\Theta(h, \Gamma_1)$ such that
\begin{itemize}
\item[(i)] $\Theta(h, \Gamma_1)$ is a region generated by the vertical
 translation $h\boldsymbol{\widehat{k}}$ of a piecewise $C^{1}$ function $\Gamma_1$ as
  given in definition \ref{Def eligible surface and generated fissure}.

\item[(ii)] $\Theta_{i}\subseteq \mathbb{R}^{3}$ are open bounded simply
connected regions for $i = 1, 2$. These stand for rock matrix regions and satisfy
%
\begin{subequations}\label{Def information of the domain}
\begin{gather}\label{Def interface boundary conditions}
%
\begin{gathered}
\partial \Theta_{1} \cap \partial\Theta(h, \Gamma_1) = \Gamma_1 ,\\
\partial \Theta_{2} \cap \partial\Theta(h, \Gamma_1)= \Gamma_1+h ,
\end{gathered}
\\
\label{Def fissures connectedness condition}
cl(\Theta_{1}) \cap cl(\Theta_{2}) = \emptyset.
\end{gather}
\end{subequations}
\end{itemize}
The sets $\Omega_{1}, \Omega_{2}$ are the rock matrix and the fissures regions
respectively; i.e.,
\begin{equation}\label{Def fissured medium region}
\Omega_{1} := \Theta_{1}\cup \Theta_{2}, \quad
\Omega_{2}:=   \Theta(h, \Gamma_1), \quad
\Omega  :=  cl(\Omega_{1}\cup\Omega_{2}).
\end{equation}
The interfaces are defined by
\begin{equation}\label{Def global interfaces}
\Gamma_1,\quad
\Gamma_2 :=  \Gamma_1+h, \quad
\Gamma := \Gamma_2\cup\Gamma_1
\end{equation}
%
The upwards normal vector to the surface $\Gamma_1$ is denoted $\boldsymbol{\widehat{n}}$; i.e.,
\begin{equation}\label{Def Standarazied upwards normal vector}
   \boldsymbol{\widehat{n}} := \frac{(-\boldsymbol{\widetilde{\nabla}} \zeta, 1)}{| (-\boldsymbol{\widetilde{\nabla}} \zeta, 1)|} .
\end{equation}
Finally, we introduce the notation $\Omega = (\Theta_{1}, \Theta_{2}, \zeta, h)$
 for this type of domains.
\end{definition}

\begin{remark}\label{Rem connectivity condition on the domains} \rm
The condition \ref{Def fissures connectedness condition} of connectivity
of the rock matrix regions only through the fissure is not required
for modeling the problem in mixed formulation as it is presented in
section \ref{Sec mixed formulation}; however it is necessary for the
 asymptotic analysis of the system. The same holds for the simple connectedness
requirement on the domains.
\end{remark}

\subsection{A local system of coordinates}\label{Sec local system of coordinates}

Several aspects of the flow through the fissure are handled more conveniently,
when the velocities are expressed in a coordinate system consistent with
the geometry of the surface that generates the crack. Let $\Gamma$ be a
surface as defined in \eqref{Def C2 surface} and $\boldsymbol{\widehat{n}}$ the upwards normal
to the surface $\Gamma$ i.e. $\boldsymbol{\widehat{n}} = \boldsymbol{\widehat{n}}(s) = \boldsymbol{\widehat{n}} (\widetilde{\mathbf{x}})$.
 Now, for each point $\widetilde{\mathbf{x}}$ we choose a local orthonormal basis in the
following way,

\begin{equation}\label{Eq local orthogonal basis} \rm
\mathcal{B}(\widetilde{\mathbf{x}}):=
\{\boldsymbol{\widehat {e}}_1(\widetilde{\mathbf{x}}), \boldsymbol{\widehat {e}}_2 (\widetilde{\mathbf{x}}), \boldsymbol{\widehat{n}}(\widetilde{\mathbf{x}})\} .
\end{equation}
Let $M = M (\widetilde{\mathbf{x}})$ be the orthogonal matrix relating the global canonical
 basis with the local one; i.e.,
\begin{subequations}\label{Def orthogonal matrix}
\begin{gather}\label{Def orthogonal matrix 1 }
M (\widetilde{\mathbf{x}})   \boldsymbol{\widehat{\imath}} = \boldsymbol{\widehat {e}}_1(\widetilde{\mathbf{x}}) , \\
\label{Def orthogonal matrix 2}
M (\widetilde{\mathbf{x}})   \boldsymbol{\widehat{\jmath}} = \boldsymbol{\widehat {e}}_2 (\widetilde{\mathbf{x}}) , \\
\label{Def orthogonal matrix 3}
M (\widetilde{\mathbf{x}})   \boldsymbol{\widehat{k}} = \boldsymbol{\widehat{n}} (\widetilde{\mathbf{x}}) .
\end{gather}
\end{subequations}
The block matrix notation for this local matrix will be
\begin{equation}\label{Eq blocks of the rotational matrix}
M (\widetilde{\mathbf{x}}) :=\begin{pmatrix}
 M^{ T,\boldsymbol{\tau}} &  M^{T,\boldsymbol{\widehat{n}}}\\
 M^{\boldsymbol{\widehat{k}},\boldsymbol{\tau}} &  M^{\boldsymbol{\widehat{k}},\boldsymbol{\widehat{n}}}
\end{pmatrix}(\widetilde{\mathbf{x}}) .
\end{equation}
Here, the index $T$ stands for the first two components in the directions
 $\boldsymbol{\widehat{\imath}},   \boldsymbol{\widehat{\jmath}}$ while the index $\tau$ stands for the tangential
velocity, orthogonal to $\boldsymbol{\widehat{n}}$. Then $\mathbf{w} =  [\mathbf{w}_{\tau}, \mathbf{w}_{\boldsymbol{\widehat{n}}}](\widetilde{\mathbf{x}})$
with the following relations
\begin{subequations}\label{Eq definition of normal and tangential velocity}
\begin{gather}\label{Eq definition of normal velocity}
w_{\boldsymbol{\widehat{n}}} := \mathbf{w}\cdot\boldsymbol{\widehat{n}} (\widetilde{\mathbf{x}}) , \\
\label{Eq definition of tangential velocity}
\mathbf{w}_{\tau} := ( \mathbf{w}\cdot \boldsymbol{\widehat {e}}_1 (\widetilde{\mathbf{x}}), \mathbf{w}\cdot\boldsymbol{\widehat {e}}_2(\widetilde{\mathbf{x}}) ) .
\end{gather}
\end{subequations}
Clearly, the relationship between velocities is given by
\begin{equation}\label{Eq local and global velocities}
\begin{aligned}
\mathbf{w} (\widetilde{\mathbf{x}}, x_{3})
&= \left\{\begin{matrix}
\widetilde{\mathbf{w}}\\
\mathbf{w}\cdot\boldsymbol{\widehat{k}}
\end{matrix}
\right\}
(\widetilde{\mathbf{x}}, x_{3})
= M (\widetilde{\mathbf{x}}) \left\{\begin{matrix}
\mathbf{w}_{\tau}\\
\mathbf{w}\cdot\boldsymbol{\widehat{n}}
\end{matrix}
\right\}
(\widetilde{\mathbf{x}}, x_{3})\\[4pt]
%
&= \begin{pmatrix}
 M^{ T,\boldsymbol{\tau}} (\widetilde{\mathbf{x}}) &  M^{T,\boldsymbol{\widehat{n}}} (\widetilde{\mathbf{x}})\\
 M^{\boldsymbol{\widehat{k}},\boldsymbol{\tau}} (\widetilde{\mathbf{x}}) &  M^{\boldsymbol{\widehat{k}},\boldsymbol{\widehat{n}}} (\widetilde{\mathbf{x}})
\end{pmatrix}
\left\{\begin{matrix}
\mathbf{w}_{\tau}\\
\mathbf{w}\cdot\boldsymbol{\widehat{n}}
\end{matrix}
\right\}
(\widetilde{\mathbf{x}}, x_{3}) .
\end{aligned}
\end{equation}

\begin{proposition}\label{Th isometry change of coordinates}
Let $h>0$, $\Upsilon$, $\Theta(h, \Upsilon)$ be as in definition
\ref{Def eligible surface and generated fissure}; let $\boldsymbol{\widehat{n}}$ be the upwards
normal to the surface $\Upsilon$ and $M$ be the matrix defined by
\eqref{Def orthogonal matrix}. Then
%
\begin{itemize}
   \item[(i)]
   The map $\mathbf{w}\mapsto M(\widetilde{\mathbf{x}}) \mathbf{w}$ is an isometry in
 $\mathbf{L}^2(\Theta(h, \Upsilon) )$. In particular if $\mathbf{w}_{\tau}, \mathbf{w}\cdot\boldsymbol{\widehat{n}}$
are defined as in \eqref{Eq definition of normal and tangential velocity}
then $\mathbf{w}\in \mathbf{L}^2(\Theta(h, \Upsilon) )$ if and only if
$\mathbf{w}_{\tau}\in L^{2}(\Theta(h, \Upsilon) )\times L^{2}(\Theta(h, \Upsilon))$
and $\mathbf{w}\cdot \boldsymbol{\widehat{n}} \in L^{2}(\Theta(h, \Upsilon))$.

   \item[(ii)]
   If $\mathbf{w}\in \mathbf{L}^2(\Theta(h, \Upsilon))$ is such that $\partial_{z} \mathbf{w} \in \mathbf{L}^2(\Theta(h, \Upsilon))$ then
%
%
\begin{equation}\label{Eq partial z velocities}
\partial_z \mathbf{w} (\widetilde{\mathbf{x}}, z)
= M (\widetilde{\mathbf{x}}) \left\{\begin{matrix}
\partial_z \mathbf{w}_{\tau}\\
(\partial_z \mathbf{w})\cdot\boldsymbol{\widehat{n}}
\end{matrix}
\right\}
(\widetilde{\mathbf{x}}, z) .
\end{equation}
\end{itemize}
\end{proposition}

\begin{proof}
(i) For $\widetilde{\mathbf{x}}$ fixed the matrix $M(\widetilde{\mathbf{x}})$ is orthogonal; i.e.,
 for arbitrary functions $\mathbf{v},\mathbf{w} \in \mathbf{L}^2(\Theta(h, \Upsilon) )$ and
 $\mathbf{x}\in \Theta(h, \Upsilon)$ holds
$\mathbf{v}(\mathbf{x})\cdot\mathbf{w}(\mathbf{x}) = M(\widetilde{\mathbf{x}}) \mathbf{v}(\mathbf{x})\cdot M(\widetilde{\mathbf{x}}) \mathbf{w}(\mathbf{x})$. Hence
   \begin{align*}%\label{Eq isometry change of coordinates}
&\int_{\Theta(h, \Upsilon) }\mathbf{v}(\mathbf{x})\cdot\mathbf{w}(\mathbf{x}) d\mathbf{x} \\
&=\int_{\Theta(h, \Upsilon) } M(\widetilde{\mathbf{x}})\mathbf{v}(\mathbf{x})\cdot M(\widetilde{\mathbf{x}})\mathbf{w}(\mathbf{x}) d\mathbf{x} \\
&=\int_{\Theta(h, \Upsilon) } \mathbf{v}_{\tau}(\mathbf{x})\cdot \mathbf{w}_{\tau}(\mathbf{x}) d\mathbf{x}
   + \int_{\Theta(h, \Upsilon) } (\mathbf{v}\cdot\boldsymbol{\widehat{n}}) (\mathbf{x})\cdot (\mathbf{w}\cdot\boldsymbol{\widehat{n}})(\mathbf{x}) d\mathbf{x}
\end{align*}
   The equality of the second line shows the necessity and sufficiency of
the tangential and normal components being square integrable in the domain
$\Theta(h, \Upsilon)$.

(ii)  It follows from a direct calculation of distributions with
$\boldsymbol {\varphi}\in [C_0^{\infty} (\Theta(h, \Upsilon) ) ]^{3}$
arbitrary and the fact that $\partial_{z} M = 0$.
\end{proof}

\subsection{The problem and its formulation}\label{Sec strong general problem}
In this section we define the problem in a rigorous way and give a variational
formulation in which it is well-posed. Let
$\Omega = (\Theta_{1}, \Theta_{2}, \zeta, h)$ be a fractured domain with a
vertical translation generated fissure as in definition
\ref{Def unidirectional translation medium}. We denote $\mathbf{v}^1,  p^1$
the velocity and pressure in the rock matrix region  $\Omega_{1}$.
In the same fashion $\mathbf{v}^2,  p^2$ denote the velocity and pressure in the
fissures region $\Omega_{2}$.

Consider the problem
\begin{subequations}\label{Pblm fissures system strong problem}
\begin{gather}\label{Eq Darcy rock matrix}
 a_1 \mathbf{u}^1 + \boldsymbol{\nabla}  p^1 + \mathbf{g} = 0\quad \text{and} \\
\label{Eq fluid mass conservation rock matrix}
\boldsymbol{\nabla}\cdot \mathbf{u}^1 = F\quad \text{in }\Omega_{1} , \\
\label{Eq drained condition rock matrix}
 p^1 = 0 \quad \text{on } \partial\Omega_{1}-\Gamma , \\
\label{Eq interface balance 1 gamma}
 p^1 -  p^2 = \alpha  \mathbf{u}^1\cdot\boldsymbol{\widehat{n}}  \boldsymbol{\mathbbm{1}}_{\Gamma_2 }
-\alpha  \mathbf{u}^1\cdot\boldsymbol{\widehat{n}}  \boldsymbol{\mathbbm{1}}_{\Gamma_1 }
\quad \text{and} \\
\label{Eq interface balance 2 gamma}
(\mathbf{u}^1 -\mathbf{u}^2 )\cdot\boldsymbol{\widehat{n}} \boldsymbol{\mathbbm{1}}_{\Gamma_1 }
-(\mathbf{u}^1 -\mathbf{u}^2 )\cdot\boldsymbol{\widehat{n}} \boldsymbol{\mathbbm{1}}_{\Gamma_2 } =
f_{\Gamma}
\quad \text{on } \Gamma , \\
\label{Eq Darcy fissures region}
 a_2 \mathbf{u}^2+\boldsymbol{\nabla} p^2 +\mathbf{g} = 0\quad \text{and} \\
\label{Eq fluid mass conservation fissures}
\boldsymbol{\nabla}\cdot \mathbf{u}^2 = F \quad \text{in } \Omega_{2} ,\\
\label{Eq non flux condition fissures}
\mathbf{u}^2 \cdot\boldsymbol{\widehat{n}} = 0 \quad \text{on}  \partial \Omega_{2}-\Gamma .
\end{gather}
\end{subequations}
The flow resistance coefficients $ a_1,  a_2$ and the fluid entry
resistance coefficient $\alpha$ are assumed to be positively bounded
from below and above, see \cite{MoralesShow2}.
In equations \eqref{Eq interface balance 1 gamma} \eqref{Eq interface balance 2 gamma}
the split of cases, $\Gamma_2$ and $\Gamma_1$, is made in order to be consistent
with the sign of the upwards normal vector $\boldsymbol{\widehat{n}}$.

\subsection{Mixed formulation of the problem}\label{Sec mixed formulation}

We start defining the spaces of velocity and pressure.
\begin{subequations}\label{Def spaces of functions}
\begin{gather}\label{Def space of velocities}
\mathbf{V}:=\{\mathbf{v}\in\mathbf{L}^2(\Omega): \boldsymbol{\nabla}\cdot\mathbf{v}^1\in \mathbf{L}^2(\Omega_{1})
 ,\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma}\in L^{2}(\Gamma)\} , \\
\label{Def space of pressures}
Q:=\{q\in L^{2}(\Omega): \boldsymbol{\nabla} q^2\in \mathbf{L}^2(\Omega_{2}) \}
\end{gather}
Endowed with their natural norms
\begin{gather}\label{Def normm space of velocities}
\|  \mathbf{v}\| _{\mathbf{V}}
:= \{\|\mathbf{v}\| _{\mathbf{L}^2(\Omega)}^{2}+\|
\boldsymbol{\nabla}\cdot\mathbf{v}^1\| _{L^{2}(\Omega_1)}^{2}
+ \|  \mathbf{v}^1\cdot\boldsymbol{\widehat{n}}\| _{L^{2}(\Gamma)}^{2}\}^{1/2},\\
\label{Def norm of pressures}
\|  q\| _{Q} := \{ \|  q\| _{L^{2}(\Omega)}^{2}
+\|  \boldsymbol{\nabla} q^2\| _{L^{2}(\Omega_{2})}^{2} \}^{1/2}
\end{gather}
\end{subequations}


\begin{remark}\label{Rem notation simplification of Gamma} \rm
In the spaces above it is understood that
\begin{equation}\label{Eq notation condensation}
   \| \mathbf{v}\cdot \boldsymbol{\widehat{n}}\|_{L^{2}(\Gamma)}^{2} =
   \| \mathbf{v}\cdot \boldsymbol{\widehat{n}}\|_{L^{2}(\Gamma_2)}^{2} +
   \| \mathbf{v}\cdot \boldsymbol{\widehat{n}}\|_{L^{2}(\Gamma_1)}^{2}\,.
   \end{equation}
\end{remark}

Consider the problem: Find $p\in Q$ and $\mathbf{u}\in \mathbf{V}$ such that
\begin{subequations}\label{Pblm weak form epsilon free system}
\begin{gather}
\label{Pblm weak form epsilon free system 1}
\begin{aligned}
&\int_{\Omega_1} a_1  \mathbf{u} \cdot \mathbf{v}
+  \int_{\Omega_2} a_2  \mathbf{u} \cdot \mathbf{v}
- \int_{\Omega_1} p  \boldsymbol{\nabla}\cdot\mathbf{v}
+ \int_{\Omega_2} \boldsymbol{\nabla} p \cdot\mathbf{v}\\
&+\alpha\int_{\Gamma}(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}})(\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}) \,dS
-\int_{\Gamma_1} p^2 (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})\,dS
+\int_{\Gamma_2} p^2 (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})\,dS\\
&= - \int_{\Omega}\mathbf{g} \cdot \mathbf{v}\,,
\end{aligned} \\
\label{Pblm weak form epsilon free system 2}
\begin{aligned}
&\int_{\Omega_1}\boldsymbol{\nabla}\cdot\mathbf{u} q - \int_{\Omega_{2}} \mathbf{u} \cdot \boldsymbol{\nabla} q
+ \int_{\Gamma_1}(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}})  q^2\,dS
- \int_{\Gamma_2}(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}})  q^2\,dS\\
&= \int_{\Omega}F q
+ \int_{\Gamma} f_{\Gamma}  q^2 \,dS\quad
\text{for all }  q\in Q,\; \mathbf{v}\in\mathbf{V}\,.
\end{aligned}
\end{gather}
\end{subequations}

\begin{remark}\label{Rem sign explanation} \rm
In the formulation above the non-symmetric interface terms are split in two pieces in order to express everything in terms of the upwards normal vector $\boldsymbol{\widehat{n}}$. In the case of the symmetric term $\int_{\Gamma}(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}})(\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})dS$ in \eqref{Pblm weak form epsilon free system 1} such split becomes unnecessary since the sign of the normal vector changes in both factors canceling each other.
\end{remark}

Define the bilinear forms $\mathcal{A}:\mathbf{V}\to\mathbf{V}'$,
$\mathcal{B}:V\to Q'$, $\mathcal{C} :Q\to Q'$ by
\begin{subequations}\label{Def operators fissures system}
\begin{gather}\label{Def velocity operator fissures system}
\mathcal{A}\mathbf{v}(\mathbf{w}):= \int_{\Omega_{1}} a_1 \mathbf{v}\cdot\mathbf{w}
+\int_{\Omega_{2}} a_2 \mathbf{v}\cdot\mathbf{w} +
\alpha \int_{\Gamma}(\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})(\mathbf{w}^1\cdot\boldsymbol{\widehat{n}}) \,dS, \\
\label{Def mixed operator fissures system}
\mathcal{B}\mathbf{v} (q) := -\int_{\Omega_{1}}\boldsymbol{\nabla}\cdot\mathbf{v} q
+\int_{\Omega_{2}}\mathbf{v}\cdot \boldsymbol{\nabla} q
-\int_{\Gamma_1} (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})   q^2 \,dS
+\int_{\Gamma_2} (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})   q^2 \,dS\,.
\end{gather}
\end{subequations}
Then, the system \eqref{Pblm weak form epsilon free system}
is a mixed formulation for the problem \eqref{Pblm fissures system strong problem}
with the abstract form
\begin{equation}\label{Pblm mixed formulation}
\begin{gathered}
\mathbf{u}\in \mathbf{V}, \;p\in Q  : \quad
\mathcal{A}\mathbf{u}+\mathcal{B} 'p = -\mathbf{g}\;\text{in }\mathbf{V} ' ,\\
-\mathcal{B} \mathbf{u} = f\quad\text{in } Q ' .
\end{gathered}
\end{equation}
For the sake of completeness we recall some well known results

\begin{theorem}\label{Th General Mixed Formulation well posedness}
   Let $\mathbf{V}, Q$ be Hilbert spaces and $\| \cdot \|_{\mathbf{V}}, \| \cdot \|_{Q}$
be their respective norms. Let $\mathcal {A}: \mathbf{V} \to \mathbf{V} '$,
$\mathcal {B}: \mathbf{V} \to Q'$ be continuous linear operators such that
   \begin{itemize}
\item[(i)]  $\mathcal {A}$ is non-negative and $\mathbf{V}$-coercive on $\ker \mathcal{B}$.

\item[(ii)] The operator $\mathcal{B}$ satisfies the inf-sup condition
   \begin{equation}\label{Ineq general inf sup condition}
         \inf_{q\in Q}   \sup_{\mathbf{v}\in \mathbf{V}} \frac{| \mathcal{ B } \mathbf{v} (q) | }{\| \mathbf{v}\|_{\mathbf{V}}   \| q\|_{Q}}  > 0 .
      \end{equation}
Then, for each $\mathbf{g}\in \mathbf{V}'$ and $f\in Q'$ there exists a unique solution
$[\mathbf{u}, p]\in \mathbf{V}\times Q$ to the problem \eqref{Pblm mixed formulation}.
Moreover, it satisfies the estimate
 \begin{equation}\label{Ineq well posedness estimate}
\| \mathbf{u}\|_{\mathbf{V}} + \| p \|_{Q}
         \leq K (\| \mathbf{g}\|_{\mathbf{V}'} + \| f\|_{Q'})\,.
      \end{equation}
   \end{itemize}
\end{theorem}

For a proof of the above theorem, see \cite{GiraultRaviartFEM} .

\begin{lemma}\label{Th control with trace and gradient}
Let $\mathcal{O}$ be an open connected bounded set in $\mathbb{R}^{ N}$
and $\mathcal{G}\subseteq \partial\mathcal{O}$ with non-null
$\mathbb{R}^{ N-1}$-Lebesgue measure, then there exists
$\kappa = \kappa(\mathcal{O}) >0$ such that
\begin{equation} \label{Eq control with trace and gradient}
\| \boldsymbol{\nabla} \eta \|_{\mathbf{L}^2(\mathcal{O})} + \|  \eta \|_{L^{2}(\mathcal{G}) }
\geq \kappa \|  \eta \|_{H^{1}(\mathcal{O})}
\end{equation}
for all $\eta\in H^{1}(\mathcal{O})$.
\end{lemma}

For a proof of the above lemma, see \cite[Proposition 5.2]{Showalter97}
or \cite[Lemma 1.2]{MoralesShow2}.

\begin{corollary}\label{Th control by gradient and trace}
   There exists a constant $\kappa >0$ such that
   \begin{equation}\label{Ineq control by gradient and trace}
      \| \boldsymbol{\nabla} q \|_{\mathbf{L}^2(\Omega_{2})} ^{2}
      + \|  q \|_{L^{2}(\Gamma)} ^{2}
      \geq\kappa
      \|  q \|_{L^{2}(\Omega_{2})} ^{2}
   \end{equation}
   for all $q\in H^{1}(\Omega_{2})$.
\end{corollary}

\begin{lemma}\label{Th inf-sup condition original}
   The operator $\mathcal {B}$ satisfies the inf-sup condition
\eqref{Ineq general inf sup condition}.
\end{lemma}

   \begin{proof}
We use the same strategy presented lemma 1.3 in \cite{MoralesShow2} with a
slight modification in the construction of the particular test function.
Fix $q\in Q$ and for $j = 1, 2$ denote $\xi_{j}$ the unique solution of
the problem
 \begin{equation}\label{Pblm Test Function}
\begin{gathered}
         -\boldsymbol{\nabla}\cdot \boldsymbol{\nabla} \xi_{j} = q^1  
\quad \text{in }  \Theta_{j}, \\
         \boldsymbol{\nabla} \xi_{j} \cdot \boldsymbol{\widehat{n}}
 =  (-1)^{j - 1}  q^2 \quad \text{on }  \Gamma_{j} ,
         \quad \boldsymbol{\nabla} \xi_{j} \cdot \boldsymbol{\widehat{n}} 
= -  q^2 \quad \text{on }
 \Gamma_{j} +  h_j \,,\\
   \xi_{j} = 0 \quad \text{on } \partial \Theta_{j} - \Gamma_{j} \,.
 \end{gathered}
      \end{equation}
Define $\mathbf{v}^1 := \boldsymbol{\nabla} \xi_{1} \boldsymbol{\mathbbm{1}}_{\Theta_{1}} 
+ \boldsymbol{\nabla} \xi_{2} \boldsymbol{\mathbbm{1}}_{\Theta_{2}}$.
Thus, $-\boldsymbol{\nabla}\cdot \mathbf{v}^1 
=   q^1  \boldsymbol{\mathbbm{1}}_{\Theta_{1}} +  q^1  \boldsymbol{\mathbbm{1}}_{\Theta_{2}}$ and
      %
      \begin{equation*}
         \mathbf{v}^1\cdot\boldsymbol{\widehat{n}} =   q^2 \boldsymbol{\mathbbm{1}}_{\Gamma_1} -   q^2 \boldsymbol{\mathbbm{1}}_{\Gamma_2}\,.
      \end{equation*}
 By the Poincar\'e inequality
$c_{1}\|  \mathbf{v}^1 \|_{H_{\rm div} (\Omega_{1})} \leq  \|  q^1 \|_{L^{2}(\Omega_{1})}
+ \|  q^2 \|_{L^{2}(\Gamma)}$. Hence, setting $\mathbf{v}^2 := \boldsymbol{\nabla}  q^2$ we have
\begin{equation}\label{Ineq inf sup inequalities chain}
\begin{aligned}
\mathcal {B}\mathbf{v} (q)
&=-\int_{\Omega_{1}}\boldsymbol{\nabla}\cdot\mathbf{v}^1  q^1 
+\int_{\Omega_{2}}\mathbf{v}^2\cdot \boldsymbol{\nabla}  q^2
-\int_{\Gamma_1} (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})   q^2 \,dS
+\int_{\Gamma_2} (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})   q^2 \,dS\\
&= \int_{\Omega_{1}}|  q^1 |^{2}
+\int_{\Omega_{2}}| \boldsymbol{\nabla}  q^2 |^{2}
+\int_{\Gamma_1} |  q^2 |^{2} \,dS
+\int_{\Gamma_2} |  q^2 |^{2} \,dS \\
&\geq \int_{\Omega_{1}}|  q^1 |^{2}
+\frac{\kappa}{2}\int_{\Omega_{2}}|   q^2 |^{2}
+ \frac{1}{2} (\int_{\Omega_{2}}|   q^2 |^{2}
+\int_{\Gamma} |  q^2 |^{2} \,dS) \\
&\geq c \| \mathbf{v}\|_{\mathbf{V}}   \| q \|_{Q}\,,
\end{aligned}
\end{equation}
for $c := \min \{c_{1}, \frac{1}{2}, \frac{\kappa}{2}\}$, which gives
the inf-sup condition of the operator $\mathcal{B}$.
   \end{proof}

\begin{theorem}\label{Th well-posedness}
Let $a_{i}(\cdot)\in L^{\infty} (\Omega)$ and
\begin{equation}\label{Def minimum coefficient coerciveness}
a^{*} := \min _{i  =  1, 2} \operatorname{ess\,inf}
 \{a_{i}(\mathbf{x}): \mathbf{x}\in \Omega_{i}\} .
\end{equation}
Then, if $a^{*}$ is positive and $0 < \alpha$, the mixed variational
formulation \eqref{Pblm mixed formulation} (or equivalently,
the system \eqref{Pblm weak form epsilon free system}) is well-posed.
\end{theorem}

\begin{proof}
Clearly $\mathcal{A}$ is non-negative and $\mathbf{V}$-coercive on $\ker \mathcal{B}$.
The operator $\mathcal {B}$ satisfies the inf-sup condition as seen in the
preceding lemma. By theorem \ref{Th General Mixed Formulation well posedness}
the result follows.
\end{proof}

\section{Scaling the problem and convergence statements}   %\label{intro}

The successful asymptotic analysis of the problem
\eqref{Pblm weak form epsilon free system} demands scaling of the heights
and resistance coefficients of the medium. We have the following definition
(see figure \ref{Fig Domains Mapping}).

\begin{definition}\label{Def epsilon parametrization fissured medium} \rm
Let $\Omega = (\Theta_{1}, \Theta_{2}, \zeta, h)$ be a fractured medium
with a vertical translation generated fissure. For $\epsilon\in (0, 1)$
we define its associated $\epsilon$-scaled fissured  system
$\{ (\Theta_{1}^{\epsilon}, \Theta_{2}^{\epsilon}, \zeta, \epsilon  h):
\epsilon>0\}$ by
\begin{subequations}\label{Eq epsilon fractured medium}
\begin{gather}\label{Eq epsilon shifts of rock matrix}
\Theta_{1}^{\epsilon} := \Theta_{1} ,\quad
\Theta_{2}^{ \epsilon} := \Theta_{2}-(1-\epsilon) h ,\\
\label{Eq epsilon fissures surfaces}
 \Gamma_1^{\epsilon} = \Gamma_1 , \quad
 \Gamma_2^{\epsilon} = \Gamma_2 -(1-\epsilon) h.
\end{gather}
\end{subequations}
The domains $\Omega_{1}^{\epsilon}, \Omega_{2}^{\epsilon}, \Omega^{\epsilon}$
and $\Gamma_1^{\epsilon}, \Gamma^{\epsilon}, \Gamma^{\epsilon}$
are defined as in \eqref{Def fissured medium region} and
\eqref{Def global interfaces} respectively.
\end{definition}

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.7\textwidth]{fig2}
\end{center}  % DomainsMappingSingleFissure.pdf
\caption{Domains mapping}\label{Fig Domains Mapping}
\end{figure}


\begin{remark} \rm
   Clearly the system $(\Theta_{1}^{\epsilon}, \Theta_{2}^{\epsilon},
\zeta, \epsilon  h)$, satisfies the conditions of definition
\ref{Def unidirectional translation medium} for every $\epsilon > 0$.
\end{remark}

\subsection{Isomorphisms of spaces and formulation}   %\label{intro}

Let $\Omega_{1}^{\epsilon}, \Omega_{2}^{\epsilon}$, $\Omega^{\epsilon}$
and $\Gamma_1^{\epsilon}, \Gamma_2^{\epsilon}$, $\Gamma^{\epsilon}$
be the domains and surfaces associated to the family
$\{ (\Theta_{1}^{\epsilon}, \Theta_{2}^{\epsilon}, \zeta, \epsilon  h):
\epsilon>0\}$ as in definition \ref{Def epsilon parametrization fissured medium}.
Define the spaces
\begin{subequations}\label{Def epsilon spaces of functions}
\begin{gather}\label{Def epsilon space of velocities}
\mathbf{V}^{\epsilon} := \{\mathbf{v}\in\mathbf{L}^2(\Omega^{ \epsilon}):
\boldsymbol{\nabla}\cdot\mathbf{v}^1\in \mathbf{L}^2(\Omega_{1}^{ \epsilon}) ,
\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma^{ \epsilon}}\in L^{2}(\Gamma^{ \epsilon}) \} , \\
\label{Def epsilon space of pressures}
Q^{\epsilon} := \{q\in L^{2}(\Omega^{\epsilon}):
\boldsymbol{\nabla} q^2\in \mathbf{L}^2(\Omega_{2}^{ \epsilon}) \} .
\end{gather}
We endow these spaces with the norms coming from the natural inner product
\begin{gather}\label{Def epsilon normm space of velocities}
\|  \mathbf{v}\| _{\mathbf{V}^{\epsilon}} :=
\{\|\mathbf{v}\| _{\mathbf{L}^2(\Omega^{\epsilon})}^{2}
+\|  \boldsymbol{\nabla}\cdot\mathbf{v}^1\| _{L^{2}(\Omega_1^{\epsilon})}^{2}
+\|  \mathbf{v}^1\cdot\boldsymbol{\widehat{n}}\| _{L^{2}(\Gamma^{\epsilon})}^{2}
\}^{1/2},\\
\label{Def epsilon norm of pressures}
\|  q\| _{Q^{\epsilon}}
:= \{\|  q\| _{L^{2}(\Omega^{\epsilon})}^{2}
+\| \boldsymbol{\nabla} q^2\| _{L^{2}(\Omega_{2}^{\epsilon})}^{2}\}^{1/2}\,.
\end{gather}
\end{subequations}
Consider the scaled problem:
Find  $p^{\epsilon}\in Q^{\epsilon}$ and $\mathbf{u}^{\epsilon}\in \mathbf{V}^{\epsilon}$ such that
\begin{subequations}\label{Pblm weak form epsilon system}
\begin{gather}
\begin{aligned}\label{Pblm weak form epsilon system 1}
&\int_{\Omega_1^{\epsilon}} a_1  \mathbf{u}^{\epsilon} \cdot \mathbf{v}  d\mathbf{y}
+  \epsilon \int_{\Omega_2^{\epsilon}} a_2  \mathbf{u}^{\epsilon} \cdot \mathbf{v} d\mathbf{y}
- \int_{\Omega_1^{\epsilon}} p^{\epsilon}  \boldsymbol{\nabla}\cdot\mathbf{v}   d\mathbf{y}
+ \int_{\Omega_2^{\epsilon}} \boldsymbol{\nabla} p^{\epsilon} \cdot\mathbf{v}  d\mathbf{y}
\\
&+\alpha\int_{\Gamma^{\epsilon}}  (\mathbf{u}^{\epsilon,1}\cdot\boldsymbol{\widehat{n}}) (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}} ) \,dS
-\int_{\Gamma_1^{\epsilon}}  p^{\epsilon,2} (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})  \,dS
+\int_{\Gamma_2^{\epsilon}}  p^{\epsilon,2}  (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}} ) \,dS\\
&= -  \int_{\Omega^{\epsilon}}  \mathbf{g}^{ \epsilon} \cdot \mathbf{v} d\mathbf{y}\,,
\end{aligned} \\
\label{Pblm weak form epsilon system 2}
\begin{aligned}
&\int_{\Omega_1^{\epsilon}}\boldsymbol{\nabla}\cdot\mathbf{u}^{\epsilon} q d\mathbf{y}
- \int_{\Omega_{2}^{\epsilon}} \mathbf{u}^{\epsilon} \cdot \boldsymbol{\nabla} q d\mathbf{y}
+ \int_{\Gamma_1^{\epsilon}}(\mathbf{u}^{\epsilon,1}\cdot\boldsymbol{\widehat{n}})  q^2\,dS
- \int_{\Gamma_2^{\epsilon}}(\mathbf{u}^{\epsilon,1}\cdot\boldsymbol{\widehat{n}})  q^2\,dS\\
&= \int_{\Omega^{\epsilon}}F^{ \epsilon} q d\mathbf{y}
+ \int_{\Gamma^{\epsilon}} f_{
\Gamma^{\epsilon} }^{ \epsilon}   q^2 \,dS \quad
\text{for all }  q\in Q^{\epsilon}, \mathbf{v}\in\mathbf{V}^{\epsilon} .
\end{aligned}
\end{gather}
\end{subequations}
Clearly, the problem \eqref{Pblm weak form epsilon system}
is well-posed since it verifies all the hypothesis of theorem \ref{Th well-posedness}.
In order to analyze the asymptotic behavior of the solution $(\mathbf{u}^{\epsilon}, p^{\epsilon})$
as $\epsilon\downarrow 0$ the geometry of the $\epsilon$-domains must be mapped
to a common domain of reference.

\subsection{The $\epsilon$-problems in a reference domain}

We introduce the change of variable (see figure \ref{Fig Domains Mapping})
 $ \boldsymbol{ \varphi } :\Omega^{ \epsilon}\to \Omega$ defined by
\begin{equation} \label{Def concise change of variable}
\begin{aligned}
\boldsymbol{\varphi }(\mathbf{y})
&:= \mathbf{y}\boldsymbol{\mathbbm{1}}_{\Theta_{1}^{\epsilon}} (\mathbf{y}) + (\widetilde{ \mathbf{y} },
y_{3}+(1-\epsilon) h  ) \boldsymbol{\mathbbm{1}}_{\Theta_j^{ \epsilon} }(\mathbf{y})\\
&\quad + ( \widetilde{ \mathbf{y} },
\frac{1}{\epsilon}(y_{3}-\zeta^{\epsilon}(\widetilde{\mathbf{y}}))
+ \zeta^{\epsilon}(\widetilde{\mathbf{y}})
+(1-\epsilon) h ) \boldsymbol{\mathbbm{1}}_{\Theta(\epsilon h, \Gamma_1^{\epsilon} ) }(\mathbf{y})
\end{aligned}
\end{equation}
Defining $(\widetilde{\mathbf{x}}, z) := \boldsymbol {\varphi } (\mathbf{y})$
the gradients are related as follows
\begin{equation}\label{Eq concise gradient structure}
   \boldsymbol{\nabla}_{\mathbf{y}} = \left\{\begin{matrix}
            \boldsymbol{\widetilde{\nabla}}_{\mathbf{x}}\\
            \partial_{z }
            \end{matrix}\right\}\boldsymbol{\mathbbm{1}}_{\Omega_1^{\epsilon} } +
   \sum_{i} \begin{bmatrix}
               I & (1 - \frac{1}{\epsilon} )   \boldsymbol{\widetilde{\nabla}}_{\mathbf{x}} \zeta (\widetilde{\mathbf{x}})\\
               0 & \frac{1}{\epsilon}
            \end{bmatrix} \left\{\begin{matrix}
            \boldsymbol{\widetilde{\nabla}}_{\mathbf{x}}\\
            \partial_{z }
            \end{matrix}\right\} \boldsymbol{\mathbbm{1}}_{\Omega_{2}^{\epsilon}  }
\end{equation}
Here, it is understood that $I$ is the identity matrix in $\mathbb{R}^{2\times 2}$.
We write $\zeta$ instead of $\zeta^{\epsilon}$ for the sake of simplicity,
recalling that both surfaces differ only by a constant of vertical translation.

\begin{theorem}\label{Th change of variable isomorphism}
Let $\boldsymbol { \varphi }: \Omega^{\epsilon}\to\Omega$ be the change
of variable defined in equation \eqref{Def concise change of variable}.
Then, the maps defined $\Phi_{1}: \mathbf{V}\to \mathbf{V}^{\epsilon}$,
 $\Phi_{2}: Q\to Q^{\epsilon}$ defined respectively by
$(\Phi_{1}\mathbf{v})(\mathbf{y}) := \mathbf{v}(\boldsymbol { \varphi } (\mathbf{y}))$ and
$(\Phi_{2} q)(\mathbf{y}) := q(\boldsymbol { \varphi } (\mathbf{y}))$ are isomorphisms.
\end{theorem}

\begin{proof}
First notice for $\mathbf{v}\in \mathbf{V}$ and $q\in Q$ the functions
$\Phi_{1}\mathbf{v}$ and $\Phi_{2} q$ are defined on $\Omega^{\epsilon}$.
Moreover, for $\ell = 1, 2$ the restriction of the change of variable
is a bijection; i.e.,
$\boldsymbol { \varphi }:\Omega^{ \epsilon}_{\ell}\to\Omega_{\ell}$ is a bijection.
 Therefore $\mathbf{v}(\cdot)\in \mathbf{L}^2(\Omega_{\ell})$ if and only if
 $\mathbf{v}(\boldsymbol { \varphi }(\cdot))\in\mathbf{L}^2(\Omega_{\ell}^{\epsilon})$
and $q (\cdot) \in L^{2}(\Omega_{\ell})$ if and only if
$q(\boldsymbol { \varphi }(\cdot))\in L^{2}(\Omega_{\ell}^{\epsilon})$.
Even more, $\boldsymbol { \varphi }:\Gamma_1^{\epsilon}\to \Gamma_1$ and
$\boldsymbol { \varphi }:\Gamma_1^{\epsilon} + \epsilon  h\to \Gamma_1 +  h$
are bijective rigid translations. Therefore, the isomorphisms
$L^{2}(\Gamma_1^{\epsilon})\simeq L^{2}(\Gamma_1)$,
$L^{2}(\Gamma_1^{\epsilon}+\epsilon  h)\simeq L^{2}(\Gamma_1+ h)$ follow.

For the isomorphism $\Phi_{1}$ take $\mathbf{v} \in \mathbf{V}$ which is equivalent to
$\mathbf{v} (\mathbf{y})\in \mathbf{L}^2(\Omega)$ and $\boldsymbol{\nabla}_{\mathbf{y}}\cdot\mathbf{v}(\mathbf{y})\in L^{2}(\Omega_{1})$.
By the previous discussion these two conditions are equivalent to
$\mathbf{v}(\boldsymbol { \varphi }(\mathbf{y})) \in \mathbf{L}^2(\Omega^{\epsilon})$ and
$\boldsymbol{\nabla}_{\mathbf{y}}\cdot\mathbf{v}(\boldsymbol { \varphi }(\mathbf{y}))
= \boldsymbol{\nabla}_{\mathbf{y}}\cdot \mathbf{v}(\mathbf{x})\in L^{2}(\Omega_{1}^{\epsilon})$.
However, equation \eqref{Def concise change of variable} yields
$\boldsymbol{\nabla}_{\mathbf{y}}\cdot\mathbf{v}(\boldsymbol{\varphi}(\mathbf{y}) ) = \boldsymbol{\nabla}_{\mathbf{y}}\cdot \mathbf{v}(\mathbf{x})
= \boldsymbol{\nabla}_{\mathbf{x}}\cdot\mathbf{v}(\mathbf{x})$ whenever $\mathbf{x}\in \Omega_{1}$; i.e.,
 $\boldsymbol{\nabla}_{\mathbf{y}}\cdot\mathbf{v}(\mathbf{y})\in L^{2}(\Omega^{ \epsilon}_{1})$ if and only
if $\boldsymbol{\nabla}_{\mathbf{x}}\cdot\mathbf{v}(\mathbf{x}) =
\boldsymbol{\nabla}_{\mathbf{x}}\cdot\mathbf{v}(\boldsymbol {\varphi} (\mathbf{y}))\in L^{2}(\Omega_{1})$ as desired.

For the map $\Phi_{2}$, the $L^{2}$-integrability condition between spaces $Q$
and $Q^{\epsilon}$ is shown using the same arguments of the first paragraph.
It remains to show the $L^{2}$-integrability condition on the gradient.
First observe that the last row in the matrix equation
\eqref{Eq concise gradient structure} implies that
$\frac{\partial}{\partial y_{3}} q(\mathbf{y})\in L^{2}(\Omega_{2}^{\epsilon})$
if and only if $\frac{\partial}{\partial z} q(\mathbf{x})\in L^{2}(\Omega_{2})$.
Second, for the derivatives in the first two directions, the equation
\eqref{Eq concise gradient structure} yields
\begin{equation*}
\frac{\partial}{\partial  y_{\ell}} q(\mathbf{y})
= \frac{\partial}{\partial x_{\ell}} q(\mathbf{x})+
(1-\frac{1}{\epsilon})\frac{\partial}{\partial  x_{\ell}}
\zeta(\mathbf{x})\;\frac{\partial}{\partial z}  q(\mathbf{x}) ,\quad \ell = 1, 2 .
\end{equation*}
Recalling the gradient of $\zeta$ is bounded, we conclude
$\frac{\partial}{\partial y_{\ell}} q(\mathbf{y})\in L^{2}(\Omega_{2}^{\epsilon})$
if an only if $\frac{\partial}{\partial x_{\ell}} q(\mathbf{x})\in L^{2}(\Omega_{2})$
for $\ell = 1, 2$. Since $\frac{\partial}{\partial  z} q(\mathbf{x})\in L^{2}(\Omega_{2})$
is immediate, the proof is complete.
\end{proof}

We are to apply the change of variable
$\boldsymbol { \varphi }: \Omega^{ \epsilon}\to \Omega$ in the problem
\eqref{Pblm weak form epsilon system}, to this end, it is more convenient
to write the system in terms of the quantities and directions which yield
estimates agreeable with the asymptotic analysis. Hence, recalling the
definition of the upwards normal vector
\eqref{Def Standarazied upwards normal vector} the following relationships hold
\begin{subequations}\label{Eq normal and tangential velocities}
\begin{gather}\label{Eq normal velocity}
|(-\widetilde{\boldsymbol{\nabla}} \zeta, 1)|\mathbf{v}\cdot\boldsymbol{\widehat{n}}
= -\mathbf{\widetilde{v}}\cdot\widetilde{\boldsymbol{\nabla}}\zeta + v_3 , \\
\label{Eq vector orthogonal to n}
(\widetilde{\mathbf{v}}, \widetilde{\mathbf{v}}\cdot\boldsymbol{\widetilde{\nabla}}\zeta)\cdot\boldsymbol{\widehat{n}} = 0
\quad \text{in}\;  \Theta( h, \Gamma_1) .
\end{gather}
\end{subequations}

Applying the change of variable \eqref{Def concise change of variable}
to  problem  \eqref{Pblm weak form epsilon system} and combining with
 relation \eqref{Eq normal velocity} we obtain the following variational statement:
Find $p^{\epsilon}\in Q$ and $\mathbf{u}^{\epsilon}\in \mathbf{V}$ such that
\begin{subequations}\label{Pblm associated weak form epsilon system}
\begin{gather}
\label{Pblm associated weak form epsilon system 1}
\begin{aligned}
&\int_{\Omega_1} a_1  \mathbf{u}^{\epsilon} \cdot \mathbf{v}
+  \epsilon^{ 2} \int_{\Omega_2} a_2  \mathbf{u}^{\epsilon} \cdot \mathbf{v}
- \int_{\Omega_1} p^{\epsilon}  \boldsymbol{\nabla}\cdot\mathbf{v}\\
&+ \int_{\Omega_{2}} \epsilon
(\widetilde{\boldsymbol{\nabla}} p^{\epsilon}+\partial_{z} p^{\epsilon} \widetilde{\boldsymbol{\nabla}} \zeta) \cdot\mathbf{\widetilde{v}}
+ \int_{\Omega_{2}}  | (-\widetilde{\boldsymbol{\nabla}} \zeta, 1 )
|  \partial_{z}p^{\epsilon}(\mathbf{v}\cdot\boldsymbol{\widehat{n}}) \\
&-\int_{\Gamma_1}p^{\epsilon,2} (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})\,dS
+\int_{\Gamma_2}p^{\epsilon,2} (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})\,dS
+\alpha\int_{\Gamma}(\mathbf{u}^{\epsilon,1}\cdot\boldsymbol{\widehat{n}})(\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}) \,dS\\
&= - \int_{\Omega_{1}}\mathbf{g}^{ \epsilon} \cdot \mathbf{v}
- \epsilon \int_{\Omega_{2}}\mathbf{g}^{ \epsilon} \cdot \mathbf{v}\,.
\end{aligned} \\
\label{Pblm associated weak form epsilon system 2}
\begin{aligned}
&\int_{\Omega_1}\boldsymbol{\nabla}\cdot\mathbf{u}^{\epsilon} q
- \int_{\Omega_{2}} \epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}
\cdot(\boldsymbol{\widetilde{\nabla}} q+\partial_{z} q \widetilde{\boldsymbol{\nabla}} \zeta)
- \int_{\Omega_{2}} | (-\widetilde{\boldsymbol{\nabla}} \zeta, 1 ) |
 (\mathbf{u}^{\epsilon,2}\cdot\boldsymbol{\widehat{n}})  \partial_{z} q \\
&+ \int_{\Gamma_1}(\mathbf{u}^{\epsilon,1}\cdot\boldsymbol{\widehat{n}})  q^2\,dS
- \int_{\Gamma_2} (\mathbf{u}^{\epsilon,1}\cdot\boldsymbol{\widehat{n}} )   q^2\,dS \\
&= \int_{\Omega_{1}}F^{ \epsilon, 1} q
+ \epsilon \int_{\Omega_{2}}F^{ \epsilon, 2} q
+ \int_{\Gamma} f_{\Gamma}^{ \epsilon}   q^2 \,dS \quad
\text{for all }  q\in Q, \mathbf{v}\in\mathbf{V}\,.
\end{aligned}
\end{gather}
\end{subequations}
Finally, by the theorem \ref{Th change of variable isomorphism}
on isomorphisms of function spaces, we conclude that the problems
\eqref{Pblm associated weak form epsilon system} and
\eqref{Pblm weak form epsilon system} are equivalent.

\subsubsection{The strong rescaled problem}\label{Sec strong epsilon problem}

The solution of the problem \eqref{Pblm associated weak form epsilon system}
is the weak solution of the system of equations
\begin{subequations}\label{Eq epsilon problem strong form}
\begin{gather}\label{Eq epsilon Darcy Omega 1}
a_1 \mathbf{u}^{\epsilon,1}+\boldsymbol{\nabla}  p^{\epsilon,1}+\mathbf{g}=0  \quad \text{and }\\
\label{Eq epsilon divergence Omega 1}
\boldsymbol{\nabla}\cdot\mathbf{u}^{\epsilon,1}=f^{ \epsilon, 1} \quad \text{in } \Omega_{1} \,.\\
\label{Eq epsilon drained rock condition rock matrix}
p^{\epsilon,1} = 0 \quad \text{on } \partial \Omega_{1}-\Gamma \,. \\
\label{Eq epsilon normal stress balance interface condition}
p^{\epsilon,1}-p^{\epsilon,2} = \alpha  \mathbf{u}^1\cdot\boldsymbol{\widehat{n}}  \boldsymbol{\mathbbm{1}}_{\Gamma_2}
-\alpha  \mathbf{u}^1\cdot\boldsymbol{\widehat{n}}  \boldsymbol{\mathbbm{1}}_{\Gamma_1}
\quad \text{and} \\
\label{Eq epsilon normal flux balance interface condition}
(\mathbf{u}^{\epsilon,1} - \mathbf{u}^{\epsilon,2})\cdot\boldsymbol{\widehat{n}}  \boldsymbol{\mathbbm{1}}_{\Gamma_1}
-(\mathbf{u}^{\epsilon,1} - \mathbf{u}^{\epsilon,2})\cdot\boldsymbol{\widehat{n}}  \boldsymbol{\mathbbm{1}}_{\Gamma_2}
= f^{\epsilon}_{\Gamma}
\quad \text{on }  \Gamma \,. \\
\label{Eq epsilon Darcy Omega 2}
\begin{bmatrix}
\epsilon a_2  \mathbf{\widetilde{u}}^{\epsilon,2}+\widetilde{\boldsymbol{\nabla}}p^{\epsilon,2}
+  (1-\frac{1}{\epsilon} )  \partial_{z} p^{\epsilon,2} \widetilde{\boldsymbol{\nabla}}\zeta
+\widetilde{\mathbf{g}}^{\epsilon} \\
\epsilon^{2} a_2   u^{\epsilon,2}_3+\partial_{z} p^{\epsilon,2}+\epsilon g^{\epsilon}_{3}
\end{bmatrix}^{T}
= \mathbf{0}  \,,\\
\label{Eq epsilon divergence Omega 2}
\boldsymbol{\nabla}\cdot (\epsilon  \mathbf{\widetilde{u}}^{\epsilon,2}, \epsilon  \mathbf{\widetilde{u}}^{\epsilon,2}\cdot \boldsymbol{\widetilde{\nabla}} \zeta
+ | (-\boldsymbol{\widetilde{\nabla}} \zeta, 1)| (\mathbf{u}^{\epsilon,2}\cdot \boldsymbol{\widehat{n}}) )
 =\epsilon F^{ \epsilon, 2}
 \quad \text{in } \Omega_{2} \,, \\
\label{Eq epsilon neumann condition on Omega 2}
\mathbf{\widetilde{u}}^{\epsilon,2}\cdot\boldsymbol{\widetilde{\nu} }=0 \quad\text{on } \partial
\Omega_{2}-\Gamma \,.
\end{gather}
\end{subequations}
As before equations \eqref{Eq epsilon normal stress balance interface condition},
\eqref{Eq epsilon normal flux balance interface condition} have the separation
of cases $\boldsymbol{\mathbbm{1}}_{\Gamma_2}, \boldsymbol{\mathbbm{1}}_{\Gamma_1}$ in order to be consistent with the
upwards normal vector $\boldsymbol{\widehat{n}}$. However, the equations
\eqref{Eq epsilon normal flux balance interface condition} and
\eqref{Eq epsilon neumann condition on Omega 2} need further clarification.
Reordering and integrating by parts the second and third summands of  equation
\eqref{Pblm associated weak form epsilon system 2}, we have
\begin{align*}
&-  \int_{\Omega_{2}}\epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}\cdot(\boldsymbol{\widetilde{\nabla}} q + \partial_{z} q
\boldsymbol{\widetilde{\nabla}} \zeta)
-  \int_{\Omega_{2}} | (-\boldsymbol{\widetilde{\nabla}} \zeta, 1)| (\mathbf{u}^{\epsilon,2}\cdot \boldsymbol{\widehat{n}}) \partial_{z} q\\
&=  - \int_{\Omega_{2}}(\epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}, \epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}\cdot \boldsymbol{\widetilde{\nabla}} \zeta
+ | (-\boldsymbol{\widetilde{\nabla}} \zeta, 1)| (\mathbf{u}^{\epsilon,2}\cdot \boldsymbol{\widehat{n}}) )\cdot\boldsymbol{\nabla} q \\
&=  \int_{\Omega_{2}}\boldsymbol{\nabla}\cdot (\epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}, \epsilon
 \mathbf{\widetilde{u}}^{\epsilon,2}\cdot \boldsymbol{\widetilde{\nabla}} \zeta + | (-\boldsymbol{\widetilde{\nabla}} \zeta, 1)| (\mathbf{u}^{\epsilon,2}\cdot \boldsymbol{\widehat{n}}) ) q\\
&\quad -  \int_{\partial \Omega_{2}}q (\epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}, \epsilon
  \mathbf{\widetilde{u}}^{\epsilon,2}\cdot \boldsymbol{\widetilde{\nabla}} \zeta
+ | (-\boldsymbol{\widetilde{\nabla}} \zeta, 1)| (\mathbf{u}^{\epsilon,2}\cdot \boldsymbol{\widehat{n}}) )\cdot\boldsymbol{\widehat{\nu}} \,dS \,,
\end{align*}
where $\boldsymbol{\widehat{\nu}}$ is the outwards unit normal vector of the boundary
$\Omega_{2}$. We focus on the boundary term
\begin{align*}
&\int_{\partial \Omega_{2} }q (\epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}, \epsilon
 \mathbf{\widetilde{u}}^{\epsilon,2}\cdot \boldsymbol{\widetilde{\nabla}} \zeta + | (-\boldsymbol{\widetilde{\nabla}} \zeta, 1)|
(\mathbf{u}^{\epsilon,2}\cdot \boldsymbol{\widehat{n}}) )\cdot\boldsymbol{\widehat{\nu}} \,dS\\
&=  \int_{\partial \Omega_{2} - (\Gamma_1 \cup \Gamma_{2}) }q
 (\epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}, \epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}\cdot \boldsymbol{\widetilde{\nabla}} \zeta
 + | (-\boldsymbol{\widetilde{\nabla}} \zeta, 1)| (\mathbf{u}^{\epsilon,2}\cdot \boldsymbol{\widehat{n}}) )\cdot\boldsymbol{\widehat{\nu}}\,dS\\
&\quad + \sum_{\ell   =   0, 1} \int_{\Gamma_{\ell} }q
 (\epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}, \epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}\cdot \boldsymbol{\widetilde{\nabla}} \zeta
+ | (-\boldsymbol{\widetilde{\nabla}} \zeta, 1)| (\mathbf{u}^{\epsilon,2}\cdot \boldsymbol{\widehat{n}}) )\cdot\boldsymbol{\widehat{\nu}}  \,dS .
\end{align*}
The equality $\boldsymbol{\widehat{\nu}}\cdot\boldsymbol{\widehat{k}} = 0$ holds on the portion of the vertical
 wall $\partial \Omega_{2} - (\Gamma_1 \cup \Gamma_{2})$; i.e.,
 the equation \eqref{Eq epsilon neumann condition on Omega 2} follows.
For the remaining pieces of the boundary recall $\boldsymbol{\widehat{n}} = \boldsymbol{\widehat{\nu}}$ on $ \Gamma_{2}$
and $\boldsymbol{\widehat{n}} = - \boldsymbol{\widehat{\nu}}$ on $ \Gamma_1$; together with the equation
\eqref{Def Standarazied upwards normal vector}, we obtain
\begin{align*}
&- \int_{\Gamma_{\ell} }q (\epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}, \epsilon
\mathbf{\widetilde{u}}^{\epsilon,2}\cdot \boldsymbol{\widetilde{\nabla}} \zeta + | (-\boldsymbol{\widetilde{\nabla}} \zeta, 1)|
(\mathbf{u}^{\epsilon,2}\cdot \boldsymbol{\widehat{n}}) )\cdot\boldsymbol{\widehat{\nu}}\,dS\\
&= (-1)^{\ell - 1}\int_{\Gamma_{\ell} }  q (\epsilon   \mathbf{\widetilde{u}}^{\epsilon,2},
\epsilon   \mathbf{\widetilde{u}}^{\epsilon,2}\cdot \boldsymbol{\widetilde{\nabla}} \zeta + | (-\boldsymbol{\widetilde{\nabla}} \zeta, 1)|
(\mathbf{u}^{\epsilon,2}\cdot \boldsymbol{\widehat{n}}) )\cdot\frac{(-\boldsymbol{\widetilde{\nabla}} \zeta, 1)}{|
(-\boldsymbol{\widetilde{\nabla}} \zeta, 1)|} \,dS \\
&= (-1)^{\ell - 1}\int_{\Gamma_{\ell} }  q (\mathbf{u}^{\epsilon,2}\cdot \boldsymbol{\widehat{n}}) \,dS
\quad\text{for }\ell = 1, 2.
\end{align*}
Combining this last identity with the interface terms in equation
\eqref{Pblm associated weak form epsilon system 2}, the strong normal
flux balance condition \eqref{Eq epsilon normal flux balance interface condition}
follows.

\subsection{A-priori estimates and convergence statements}
\label{Sec a priori estimates and convergence}

To obtain a-priori estimates on the norm of the solutions the following
hypothesis are assumed:
\begin{subequations}\label{Eq boundedness forcing terms}
\begin{gather}\label{Eq boundedness forcing divergence terms}
\| F^{ \epsilon}\|_{L^{2}(\Omega)} \text{ is bounded and }
F^{1,\epsilon}\overset{w}\rightharpoonup F^{ 1}\text{ in }L^{2}(\Omega_{1}) ,\\
\label{Eq boundedness forcing gravity terms}
\mathbf{g}^{\epsilon}\overset{w}\rightharpoonup \mathbf{g}\text{ in } \mathbf{L}^2(\Omega_{1}) ,\quad
\mathbf{g}^{ 2,\epsilon}(\widetilde{\mathbf{x}}, \epsilon z)\overset{w}\rightharpoonup
\mathbf{g}(\widetilde{\mathbf{x}})\text{ in } \mathbf{L}^2(\Omega_{2}) ,\\
\label{Eq boundedness interface forcing terms}
\text{and } f^{ \epsilon}_{
\Gamma}\overset{w}\rightharpoonup f_{\Gamma}
\text{ in }L^{2}(\Gamma) .
\end{gather}
\end{subequations}
Now test equation \eqref{Pblm associated weak form epsilon system 1}
with $\mathbf{u}^{\epsilon}$ and equation \eqref{Pblm associated weak form epsilon system 2}
with $p^{\epsilon}$, add them together and obtain
\begin{equation} \label{Eq first estimate}
\begin{aligned}
&a^{*} (\| \mathbf{u}^{\epsilon,1}\|_{0, \Omega_1}^{ 2}
+\|\epsilon \mathbf{u}^{\epsilon,2}\|_{0, \Omega_2}^{ 2})
+\alpha \|\mathbf{u}^{\epsilon,1}\cdot\boldsymbol{\widehat{n}}\|_{L^{2}(\Gamma)}^{ 2}\\
&=\int_{\Omega_{1}}F^{ \epsilon, 1} p^{\epsilon}
+\epsilon\int_{\Omega_{2}}F^{ \epsilon, 2} p^{\epsilon}
+\int_{\Gamma}f^{ \epsilon}_{\Gamma} p^{\epsilon,2}\,dS
-\int_{\Omega_1}\mathbf{g}^{ 1}\cdot\mathbf{u}^{\epsilon}
-\int_{\Omega_2}\mathbf{g}^{ 2}\cdot \epsilon \mathbf{u}^{\epsilon}\\
&\leq C (\| F^{ \epsilon}\|_{0, \Omega}
+\| f^{ \epsilon}_{\Gamma}\|_{0, \Gamma})\|  p^{\epsilon}\|_{Q}
+\| \mathbf{g}^{ \epsilon}\|_{0, \Omega}(\| \mathbf{u}^{\epsilon,1} \|_{0, \Omega_1}
+\| \epsilon \mathbf{u}^{\epsilon,2} \|_{0, \Omega_{2}}) .
\end{aligned}
\end{equation}
Here, the constant $C>0$ is independent from $\epsilon>0$.
Next, the term $\| p^{\epsilon} \|_{Q}$ must be bounded in terms of the
flux $\mathbf{u}^{\epsilon,1} \boldsymbol{\mathbbm{1}}_{\Omega_{1}} + \epsilon   \mathbf{u}^{\epsilon,2}\boldsymbol{\mathbbm{1}}_{\Omega_{2}}$
and the forcing terms. By the third component of the vector
 equation \eqref{Eq epsilon Darcy Omega 2} we have
\begin{subequations}\label{Eq boundedness pressure gradient omega 2}
\begin{equation}\label{Eq boundedness normal pressure gradient omega 2}
\| \frac{1}{\epsilon} \partial_{z}  p^{\epsilon,2} \|_{0, \Omega_{2}}
\leq \epsilon \|  a_{2}\|_{L^{ \infty}(\Omega_{2})}
\|  u^{\epsilon,2}_3 \|_{0, \Omega_{2}}+\| g^{\epsilon}_{3} \|_{0, \Omega_{2}} .
\end{equation}
Combined with the first two components of the vector equation
\eqref{Eq epsilon Darcy Omega 2} yields
\begin{equation}\label{Eq boundedness tangential pressure gradient omega 2}
\| \boldsymbol{\widetilde{\nabla}} p^{\epsilon,2} \|_{0, \Omega_{2}} \leq C (
\| a_{2}\|_{L^{ \infty}(\Omega_{2})}
\| \epsilon  \mathbf{u}^{\epsilon,2} \|_{0, \Omega_{2}}
+\| \mathbf{g}^{ \epsilon} \|_{0, \Omega_{2}}) ,
\end{equation}
\end{subequations}
for an adequate constant $C>0$. Thus
\begin{equation}\label{Eq boundedness pressure full gradient omega 2}
\| \boldsymbol{\nabla} p^{\epsilon,2} \|_{0, \Omega_{2}} \leq
C (\| a_{2}\|_{L^{\infty}(\Omega_{2})}
\| \epsilon \mathbf{u}^{\epsilon,2} \|_{0, \Omega_{2}}
+\| \mathbf{g} \|_{0, \Omega_{2}}).
\end{equation}
With $C>0$ a constant independent from $\epsilon>0$.
Additionally, the equation \eqref{Eq epsilon Darcy Omega 1} yields
\begin{equation}\label{Eq boundedness pressure full gradient omega 1}
\| \boldsymbol{\nabla} p^{\epsilon,2} \|_{0, \Omega_{1}} \leq  \| a_{1}\|_{L^{\infty}(\Omega_{2})}
\|  \mathbf{u}^{\epsilon,2} \|_{0, \Omega_{1}} +\| \mathbf{g} \|_{0, \Omega_{1}} .
\end{equation}
The boundary condition \eqref{Eq epsilon drained rock condition rock matrix}
together with Poincar\'e inequality give the control
$\|  p^{\epsilon,1} \|_{1, \Omega_{1}} \leq C \| \boldsymbol{\nabla} p^{\epsilon} \|_{0, \Omega_{1}}$.
On the other hand, the inequality \eqref{Ineq control by gradient and trace}
implies $\| p^{\epsilon} \|_{1, \Omega_{2}} \leq C (\| p^{\epsilon} \|_{0, \Gamma}
+ \| \boldsymbol{\nabla} p^{\epsilon} \|_{0, \Omega_{2}}) $; combined with the normal stress
balance conditions \eqref{Eq epsilon normal stress balance interface condition}
we conclude:
\begin{equation}\label{Ineq presssure norm controled by pressure gradient on the solution}
\| p^{\epsilon} \|_{Q} \leq \| p^{\epsilon} \|_{1, \Omega} \leq C\|  \boldsymbol{\nabla} p^{\epsilon} \|_{0, \Omega} .
\end{equation}
And $C>0$ is independent from $\epsilon>0$. Finally, a combination of
inequalities
\eqref{Ineq presssure norm controled by pressure gradient on the solution},
\eqref{Eq boundedness pressure full gradient omega 1} and
\eqref{Eq boundedness pressure full gradient omega 2} imply that the left
hand side of inequality \eqref{Eq first estimate} is bounded.
From the observations above we conclude that the following sequences are bounded
\begin{subequations}\label{Eq boundedness of epsilon solutions}
\begin{gather}\label{Eq boundedness epsilon velocities}
\| \mathbf{u}^{\epsilon,1} \|_{0, \Omega_{1}}  ,\quad
\| \epsilon \mathbf{u}^{\epsilon,2} \|_{0, \Omega_{2}} ,\quad
\sqrt{\alpha} \| \mathbf{u}^{\epsilon,1}\cdot\boldsymbol{\widehat{n}} \|_{L^{2}(\Gamma)},\\
\label{Eq boundedness of epsilon pressures}
\| p^{\epsilon,1} \|_{H^{1}(\Omega_{1})}  ,\quad
\| p^{\epsilon,2} \|_{H^{1}(\Omega_{2})} ,\quad
\| \frac{1}{\epsilon} \partial_{z} p^{\epsilon} \|_{0, \Omega_{2}} ,\quad
\| \boldsymbol{\nabla}\cdot \mathbf{u}^{\epsilon,1} \|_{L^{2}(\Omega_{1})} .
\end{gather}
\end{subequations}


\begin{remark}\label{Rem Divergence Behavior} \rm
The change of variable $\boldsymbol { \varphi }$ modifies the structure
of the divergence on the domain $ \Omega_{2}$, therefore it can only be
claimed that the linear combination $\epsilon \widetilde{\boldsymbol{\nabla}}\cdot\mathbf{\widetilde{u}}^{\epsilon,2}
+\epsilon  (1-\frac{1}{\epsilon} )\partial_{z} (\widetilde{\boldsymbol{\nabla}}\zeta\cdot\mathbf{\widetilde{u}}^{\epsilon,2} )
+ \partial_{z} u^{\epsilon,2}_3$ is bounded in $L^{2}(\Omega_{2})$.
\end{remark}

\subsection{Weak limits}\label{Sec weak convergence}

In the previous section bounds independent from $\epsilon>0$ are obtained
for $[\mathbf{u}^{\epsilon,1}, \epsilon \mathbf{u}^{\epsilon,2}]\in \mathbf{V}$ and
$p^{\epsilon} = [p^{\epsilon,1}, p^{\epsilon,2}]\in H^{1}(\Omega_{1})\times H^{1}(\Omega_{2})$,
consequently in $Q$. Then, there must exist $\mathbf{u}\in \mathbf{V}$, $p\in Q$,
$\eta\in L^{2}(\Omega_{2})$ and a subsequence, from now on denoted the same,
such that
\begin{subequations}\label{Eq convergence of epsilon solutions}
\begin{gather}\label{Eq convergence of epsilon pressures}
p^{\epsilon} \overset{w}\rightharpoonup p\quad \text{in $Q$ and strongly in }L^{2}(\Omega),\\
\label{Eq convergence of epsilon velocities 1}
\mathbf{u}^{\epsilon,1} \overset{w}\rightharpoonup \mathbf{u}^1 \text{ in } \mathbf{L}^2(\Omega_{1}) \quad
\text{and} \quad
\boldsymbol{\nabla}\cdot\mathbf{u}^{\epsilon,1} \overset{w}\rightharpoonup \boldsymbol{\nabla}\cdot\mathbf{u}^1 \text{ in } L^{2}(\Omega_{1}) ,\\
\label{Eq convergence of epsilon velocities on interface}
\sqrt{\alpha} \mathbf{u}^{\epsilon,1}\cdot\boldsymbol{\widehat{n}} \overset{w}\rightharpoonup \sqrt{\alpha} \mathbf{u}^1\cdot\boldsymbol{\widehat{n}}\quad
 \text{in } L^{2}(\Gamma) , \\
\label{Eq convergence of epsilon velocities 2}
\epsilon  \mathbf{u}^{\epsilon,1} \overset{w}\rightharpoonup \mathbf{u}^2\quad \text{in }\mathbf{L}^2(\Omega_{2}),\\
\label{Eq convergence of higher order pressure}
\frac{1}{\epsilon} \partial_{z}p^{\epsilon,2} \overset{w}\rightharpoonup \eta
 \text{ in } L^{2}(\Omega_{2}) \quad \text {and} \quad
\partial_{z}p^{\epsilon,2} \to 0 \text{ strongly in } L^{2}(\Omega_{2}) .
\end{gather}
\end{subequations}
Choose $\phi\in C_0^{\infty}(\Omega_{2})$ arbitrary, test the
equation \eqref{Pblm associated weak form epsilon system 2} with
$q := \epsilon \phi$ and let $\epsilon\downarrow 0$.
Recalling \eqref{Eq convergence of epsilon velocities 2} this gives
\begin{align*}
0 &= \lim_{\epsilon\downarrow 0} \int_{\Omega_{2}} |(-\boldsymbol{\widetilde{\nabla}} \zeta, 1)
|   (\epsilon \mathbf{u}^{\epsilon,2}\cdot\boldsymbol{\widehat{n}}) \partial_{z}\phi \\
&= \int_{\Omega_{2}} |(-\boldsymbol{\widetilde{\nabla}} \zeta, 1) |
(\mathbf{u}^2\cdot\boldsymbol{\widehat{n}}) \partial_{z}\phi\\
&= -\langle \partial_{z} |(-\boldsymbol{\widetilde{\nabla}} \zeta, 1) |
  (\mathbf{u}^2\cdot\boldsymbol{\widehat{n}}), \phi\rangle _{D '(\Omega_{2}), D(\Omega_{2})} .
\end{align*}
Since $(-\boldsymbol{\widetilde{\nabla}} \zeta, 1)$ does not depend on the vertical variable
$z$ and it is the non-zero vector almost everywhere, we conclude that
$\partial_{z}(\mathbf{u}^2\cdot\boldsymbol{\widehat{n}}) = 0$ i.e. the component of the velocity
normal to the surface $\Gamma_1$ is independent from $z$ in $\Omega_{2}$.
Now choose $q\in Q$ arbitrary, test
\eqref{Pblm associated weak form epsilon system 2} with $\epsilon  q$
and let $\epsilon \downarrow 0$ to get
\begin{align*}
0 &= \int_{\Omega_{2}} |(-\boldsymbol{\widetilde{\nabla}} \zeta, 1) |
  (\mathbf{u}^2\cdot\boldsymbol{\widehat{n}}) \partial_{z} q \; d\mathbf{x} \\
&= \int_{G}\int_{\zeta(\widetilde{\mathbf{x}})}^{ \zeta(\widetilde{\mathbf{x}})+ h}
|(-\boldsymbol{\widetilde{\nabla}} \zeta, 1)|   (\mathbf{u}^2\cdot\boldsymbol{\widehat{n}}) \partial_{z} q\; dz d\widetilde{\mathbf{x}}\\
&=   \int_{G} |(-\boldsymbol{\widetilde{\nabla}} \zeta, 1) |(\mathbf{u}^2\cdot\boldsymbol{\widehat{n}})
 [ q (\widetilde{\mathbf{x}}, \zeta(\widetilde{\mathbf{x}})+ h) - q(\widetilde{\mathbf{x}}, \zeta(\widetilde{\mathbf{x}}))]\, d\widetilde{\mathbf{x}} .
\end{align*}
The above holds for all $q\in Q$, in particular choosing
$q(\widetilde{\mathbf{x}}, \zeta(\widetilde{\mathbf{x}})) = \phi(\widetilde{\mathbf{x}})$ for $\phi\in C_0^{\infty}(G)$
 arbitrary and $q(\widetilde{\mathbf{x}}, \zeta(\widetilde{\mathbf{x}})+ h) = 0$ the statement above
transforms in
\begin{equation*}
\int_{G} |(-\boldsymbol{\widetilde{\nabla}} \zeta, 1) |
(\mathbf{u}^2\cdot\boldsymbol{\widehat{n}})(\widetilde{\mathbf{x}}, \zeta(\widetilde{\mathbf{x}}))  \phi (\widetilde{\mathbf{x}}, \zeta(\widetilde{\mathbf{x}}))
\, d\widetilde{\mathbf{x}} \quad
\forall \phi\in C_0^{\infty}(G) .
\end{equation*}
Therefore, $|(-\widetilde{\boldsymbol{\nabla}} \zeta, 1)|
(\mathbf{u}^2\cdot\boldsymbol{\widehat{n}})$ must be null and since $|(-\widetilde{\boldsymbol{\nabla}} \zeta, 1 ) |$
is non-zero almost everywhere we conclude
\begin{equation}\label{Eq null normal velocity}
\mathbf{u}^2\cdot\boldsymbol{\widehat{n}} = 0 \; \text{in} \; \Omega_{2} .
\end{equation}
This implies that the Cartesian coordinates of $\mathbf{u}^2$ satisfy the relation
\begin{equation}\label{Eq structure of limit velocity}
\mathbf{u}^2 = \left\{
\begin{matrix}
\mathbf{\widetilde{u}}^2 \\
 u_3^2
\end{matrix}\right\}=
\left\{
\begin{matrix}
\mathbf{\widetilde{u}}^2 \\
\mathbf{\widetilde{u}}^2\cdot\widetilde{\boldsymbol{\nabla}}\zeta
\end{matrix}\right\}\quad \text{in } \Omega_{2}  .
\end{equation}
Now take a function $\mathbf{v}_{\tau}^2\in(C_0^{\infty}(\Omega_{2})^{2}$.
 Recalling \eqref{Eq local and global velocities} define
$\mathbf{\widetilde{v}}:=  M^{ T,\boldsymbol{\tau}} \mathbf{v}_{\tau}^2$ and $v_3:=  M^{\boldsymbol{\widehat{k}},\boldsymbol{\tau}} \mathbf{v}_{\tau}^2$.
Then, the function $ \mathbf{v}^2:= \frac{1}{\epsilon} (\mathbf{\widetilde{v}}, v_3)$
has the structure \eqref{Eq structure of limit velocity} or equivalently
 $\mathbf{v}^2\cdot\boldsymbol{\widehat{n}} = 0$ inside $\Omega_{2}$. Define $\mathbf{v}$ as the trivial extension
of $\mathbf{v}^2$ to the whole domain $\Omega$, therefore $\mathbf{v}\in \mathbf{V}$.
Test \eqref{Pblm associated weak form epsilon system 1} with $\mathbf{v}$ and let
 $\epsilon \downarrow 0$, this gives
\[
\int_{\Omega_{2}}a_2 (\mathbf{x}) \mathbf{u}^2\cdot(\mathbf{\widetilde{v}}, v_3)
+\int_{\Omega_{2}}\widetilde{\boldsymbol{\nabla}}  p^2
\cdot\mathbf{\widetilde{v}}
+\int_{\Omega_{2}}\mathbf{g}\cdot(\mathbf{\widetilde{v}}, v_3)  = 0 .
\]
Consequently,
\[
\int_{\Omega_{2}}a_2 (\mathbf{x}) \mathbf{u}_{\tau}^2\cdot\mathbf{v}_{\tau}
+\int_{\Omega_{2}}\widetilde{\boldsymbol{\nabla}}  p^2
\cdot M^{ T,\boldsymbol{\tau}}\mathbf{v}_{\tau}
+\int_{\Omega_{2}}\mathbf{g}_{\tau}\cdot \mathbf{v}_{\tau} = 0 .
\]
The equation above holds for all $\mathbf{v}_{\tau} \in (C_0^{\infty}(\Omega_{2})^{2}$
and due to the isomorphism of proposition \ref{Th isometry change of coordinates}
we conclude
\begin{equation}\label{Eq lower dimensional Darcy-type}
a_2 (\mathbf{x}) \mathbf{u}_{\tau}^2+ ( M^{ T,\boldsymbol{\tau}})'  \widetilde{\boldsymbol{\nabla}}  p^2
+\mathbf{g}_{\tau} = 0  \quad \text{in } \Omega_{2} .
\end{equation}
Equation \eqref{Eq convergence of higher order pressure} implies that
$ p^2$ does not depend on the variable $z$ on $\Omega_{2}$; i.e.,
 $ p^2 =  p^2(\widetilde{\mathbf{x}})$. Therefore assuming
\begin{equation}\label{Eq dependence of limit data}
a_2 = a_{2}(\widetilde{\mathbf{x}}) ,\quad
\widetilde{\mathbf{g}} = \widetilde{\mathbf{g}}(\widetilde{\mathbf{x}})  \quad\text{in } \Omega_{2}
\end{equation}
the equation \eqref{Eq lower dimensional Darcy-type} gives
$\mathbf{u}_{\tau}^2 = \mathbf{u}_{\tau}^2(\widetilde{\mathbf{x}})$ i.e. $\mathbf{u}_{\tau}^2$ is independent from $z$ in
$\Omega_{2}$. Together with the fact $\mathbf{u}^2\cdot\boldsymbol{\widehat{n}} = 0$ in $\Omega_{2}$
we conclude that the whole vector velocity $\mathbf{u}^2$ is independent from
$z$ in $\Omega_{2}$.

\begin{remark}\label{Rem lower dimensional Darcy's law} \rm
   Observe that by the assumptions for the data
\eqref{Eq dependence of limit data} the equation
\eqref{Eq lower dimensional Darcy-type} is independent from $z$,
becoming a lower-dimensional Darcy-type constitutive law on the stream
lines parallel to $\zeta$.
\end{remark}

\section{The limit problem}\label{Sec limit problem}

Define the subspaces
\begin{subequations}\label{Def limiti space of functions}
\begin{gather}\label{Def limit space of velocities}
\mathbf{V}_0 := \{\mathbf{v}\in \mathbf{V}: \partial_{z}\mathbf{v}^2 = 0 ,\; \mathbf{v}^2\cdot\boldsymbol{\widehat{n}} = 0 \text{ in }
\Omega_{2} \}; \\
\label{Def limit space of pressures}
Q_0 := \{q\in Q: \partial_{z} q = 0\text{ in } \Omega_{2}\} .
\end{gather}
\end{subequations}
From the structure of the space, if $\mathbf{v} = [\mathbf{v}^1, \mathbf{v}^2]\in \mathbf{V}_0$
 then the function $[\mathbf{v}^1, \frac{1}{\epsilon} \mathbf{v}^2]$ is also in $\mathbf{V}_0$.
Hence, we use $[\mathbf{v}^1, \frac{1}{\epsilon} \mathbf{v}^2]$ to test
 \eqref{Pblm associated weak form epsilon system 1} and $q\in Q_0$ to
test \eqref{Pblm associated weak form epsilon system 2}, then we let
$\epsilon\downarrow 0$ and conclude that the limits
$[\mathbf{u}^{\epsilon,1}, \epsilon \mathbf{u}^{\epsilon,2}]\to \mathbf{u}$ and $p^{\epsilon}\to p$ are a solution of
 the \emph{limit problem}:
Find $p\in Q_0$ and $\mathbf{u} \in \mathbf{V}_0$ such that
\begin{subequations}\label{Pblm limit system}
\begin{gather}\label{Pblm limit system 1}
\begin{aligned}
&\int_{\Omega_1} a_1  \mathbf{u} \cdot \mathbf{v}
- \int_{\Omega_1} p  \boldsymbol{\nabla}\cdot\mathbf{v}
+  \int_{\Omega_2} a_2  \mathbf{u}_{\tau}^2 \cdot \mathbf{v}_{\tau}^2
+ \int_{\Omega_{2}}
\widetilde{\boldsymbol{\nabla}} p^{\epsilon} \cdot\mathbf{\widetilde{v}}  \\
&-\int_{\Gamma_1}   p^2 ( \mathbf{v}^1\cdot\boldsymbol{\widehat{n}} )  \,dS
+\int_{\Gamma_2}   p^2 ( \mathbf{v}^1\cdot\boldsymbol{\widehat{n}} )  \,dS
+\alpha  \int_{\Gamma} ( \mathbf{u}^1\cdot\boldsymbol{\widehat{n}}) (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}} ) \,dS \\
&= - \int_{\Omega_{1}}  \mathbf{g} \cdot \mathbf{v}
-  \int_{\Omega_{2}}  \mathbf{g}_{\tau} \cdot \mathbf{v}_{\tau}^2\,,
\end{aligned} \\
\label{Pblm limit system 2}
\begin{aligned}
&\int_{\Omega_1}\boldsymbol{\nabla}\cdot\mathbf{u}\; q
- \int_{\Omega_{2}} \mathbf{\widetilde{u}}^2 \cdot \boldsymbol{\widetilde{\nabla}} q
+ \int_{\Gamma_1}(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}})  q^2\,dS
- \int_{\Gamma_2}(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}})  q^2\,dS \\
&= \int_{\Omega_{1}}F^{  1} q
+ \int_{\Gamma} f_{\Gamma}   q^2 \,dS \quad
\text{for all }  q\in Q_0, \mathbf{v}\in\mathbf{V}_0 \,.
\end{aligned}
\end{gather}
\end{subequations}

\subsection{Well-Posedness of the limit problem}
\label{Sec well posedness limit problem}

Problem \eqref{Pblm limit system} is a mixed formulation of the type
\eqref{Pblm mixed formulation} with the operators $\mathcal{A}^{0}:\mathbf{V}_0\to\mathbf{V}_0'$
and $\mathcal{B }^{0}: \mathbf{V}_0 \to Q_0 '$ defined by
\begin{subequations}\label{Def limit bilinear forms}
\begin{gather}\label{Def limit bilinear form velocities}
\mathcal{A}^{0}\mathbf{v}(\mathbf{w}) := \int_{\Omega_{1}} a_1 \mathbf{v}\cdot\mathbf{w}
+\int_{\Omega_{2}} a_2 \mathbf{v}_{\tau}^2\cdot\mathbf{w}_{\tau}^2  +
\alpha \int_{\Gamma}(\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})(\mathbf{w}^1\cdot\boldsymbol{\widehat{n}}) \,dS,\\
\label{Def limit bilinear form mixed}
\mathcal{B}^{ 0}\mathbf{v}(q) := -\int_{\Omega_{1}}\boldsymbol{\nabla}\cdot\mathbf{v} \; q
+\int_{\Omega_{2}}\mathbf{\widetilde{v}}\cdot \boldsymbol{\widetilde{\nabla}} q
-\int_{\Gamma_1} (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})   q^2
+\int_{\Gamma_2} (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})   q^2\,.
\end{gather}
\end{subequations}

\begin{theorem}\label{Th limit inf-sup condition}
The operator $\mathcal{B}^{ 0}$ satisfies the inf-sup condition.
\end{theorem}

\begin{proof}
The proof has the same structure as lemma \ref{Th inf-sup condition original},
there is only one detail to be examined in the construction of the test functions.
Fix $q = [ q^1 ,  q^2]\in Q_0$, construct $\mathbf{v}^1$ in the same way it is
done in problem \eqref{Pblm Test Function}. On the other hand since
$ q^2 \in H^{ 1}(\Omega_{2})$ and $\partial_z   q^2 = 0$, define
\begin{equation*}
\mathbf{v}^2:= (\widetilde{\boldsymbol{\nabla}}  q^2, \mathbf{\widetilde{v}}\cdot\widetilde{\boldsymbol{\nabla}} \zeta )
\quad \text{in } \Omega_{2}.
\end{equation*}
Then $\mathbf{v}^2\cdot\boldsymbol{\widehat{n}} = 0$ and $\partial_{z} \mathbf{v}^2 = 0$ in $\Omega_{2}$; i.e.,
 $\mathbf{v}^2\in \mathbf{V}_0$ and $\| \mathbf{v}^2\| _{0,\Omega_2} \leq C \|  q^2\| _{1,\Omega_2}$
as desired. Repeating the inequalities presented in
\eqref{Ineq inf sup inequalities chain} the proof follows.
\end{proof}


Since the inf-sup condition holds, the theorem \ref{Th well-posedness}
applies to the operators \eqref{Def limit bilinear forms} on the spaces
$\mathbf{V}_0, Q_0$ and the limit problem \eqref{Pblm limit system} is well-posed.
By the uniqueness of the solution of the limit problem, it follows that the
original sequence converges weakly to the limit 
$\mathbf{u}\in \mathbf{V}_0,  p\in Q_0$.

\subsection{The strong form}\label{Sec strong limit problem}

To describe the \emph{strong limit problem} corresponding to
\eqref{Pblm limit system}, two properties have to be exploited.
 First, the structure $\mathbf{v}\cdot\boldsymbol{\widehat{n}} = 0$ in $\Omega_{2}$ for all
$\mathbf{v}\in \mathbf{V}_0$, this implies that $\mathbf{\widetilde{v}} =  M^{ T,\boldsymbol{\tau}} \mathbf{v}_{\tau}^2$, with $ M^{ T,\boldsymbol{\tau}}$
the matrix defined in \eqref{Eq local and global velocities}.
Second, the independence of the velocities and pressures with respect to
$z$ in $\Omega_{2}$. This last property allows to write the integrals over
$\Omega_{2}$ as surface integrals on $\Gamma_1$.
Hence, the system \eqref{Pblm limit system} transforms into:
Find $p\in Q_0$ and $\mathbf{u}\in \mathbf{V}_0$ such that
\begin{subequations}\label{Pblm lower dimensional variational statement}
\begin{gather}\label{Pblm lower dimensional variational statement 1}
\begin{aligned}
&\int_{\Omega_1} a_1  \mathbf{u} \cdot \mathbf{v}
- \int_{\Omega_1} p  \boldsymbol{\nabla}\cdot\mathbf{v}
+  h  \int_{\Gamma_1}  (\boldsymbol{\widehat{n}}\cdot\boldsymbol{\widehat{k}} )(a_2  \mathbf{u}_{\tau}^2 + ( M^{ T,\boldsymbol{\tau}})'
\widetilde{\boldsymbol{\nabla}} p^{\epsilon} +\mathbf{g}_{\tau})\cdot \mathbf{v}_{\tau}^2  \,dS\\
&-\int_{\Gamma_1} p^2 (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})\,dS
+\int_{\Gamma_2} p^2 (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}})\,dS
+\alpha\int_{\Gamma}(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}})(\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}) \,dS\\
&= - \int_{\Omega_{1}}\mathbf{g} \cdot \mathbf{v} \,,
\end{aligned}\\
\label{Pblm lower dimensional variational statement 2}
\begin{aligned}
&\int_{\Omega_1}\boldsymbol{\nabla}\cdot\mathbf{u}  q
-   h\int_{\Gamma_1} (\boldsymbol{\widehat{n}}\cdot\boldsymbol{\widehat{k}} )
 M^{ T,\boldsymbol{\tau}}\mathbf{u}_{\tau}^2 \cdot \boldsymbol{\widetilde{\nabla}} 
 q^2  \,dS
+ \int_{\Gamma_1}(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}})  q^2\,dS
- \int_{\Gamma_2}(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}})  q^2\,dS \\
&= \int_{\Omega_{1}}F^{  1} q
+ \int_{\Gamma} f_{\Gamma}   q^2\, dS \quad
\text{for all } q\in Q_0 ,\; \mathbf{v}\in \mathbf{V}_0\,.
\end{aligned}
\end{gather}
\end{subequations}
Integrating by parts the statement above, we obtain the strong problem
 with piecewise $C^{1}$ surface interfaces:
\begin{subequations}\label{Pblm limit problem strong form}
\begin{gather}\label{Pblm strong limit Darcy Omega 1}
a_1 \mathbf{u}+\boldsymbol{\nabla}  p^1+\mathbf{g}^{1}=0, \\
\label{Pblm strong limit divergence Omega 1}
\boldsymbol{\nabla}\cdot\mathbf{u}=F^{ 1} \quad \text{in } \Omega_{1}\,, \\
 p^1 = 0 \quad\text{on } \partial \Omega_{1}-\Gamma \,, \\
\mathbf{u}^2\cdot\boldsymbol{\widehat{n}}=0 ,
\quad \partial_z  p^2=0 \,,\\
\label{Pblm strong limit Darcy Omega 2}
[ a_2 (s) \mathbf{u}^2_{\tau}+( M^{ T,\boldsymbol{\tau}} )' \widetilde{\boldsymbol{\nabla}} p^2 +\mathbf{g}_{\tau}^{ 2} (s)
] \boldsymbol{\mathbbm{1}}_{\Gamma_1}= 0 \,,\\
\label{Pblm strong limit divergence Omega 2}
[ (\mathbf{u}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_{2}}-\mathbf{u}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1}) ]\boldsymbol{\mathbbm{1}}_{\Gamma_1}
+  h (\boldsymbol{\widehat{n}}\cdot\boldsymbol{\widehat{k}} ) \widetilde{\boldsymbol{\nabla}}\cdot( M^{ T,\boldsymbol{\tau}}\mathbf{u}^2_{\tau})
 \boldsymbol{\mathbbm{1}}_{\Gamma_1}
= f_{\Gamma} \quad\text{in } \Gamma \,,\\
\label{Pblm strong normal stress balance}
 p^1 -  p^2 = \alpha   \mathbf{u}^1\cdot\boldsymbol{\widehat{n}}   \boldsymbol{\mathbbm{1}}_{\Gamma_2}
-\alpha   \mathbf{u}^1\cdot\boldsymbol{\widehat{n}}    \boldsymbol{\mathbbm{1}}_{\Gamma_1} \,,\\
\label{Pblm strong limit flux boundary condition}
\mathbf{u}^2\cdot\boldsymbol{\widehat{\nu}} =  0 \quad \text{on } \partial G \,.
\end{gather}
\end{subequations}
The statement of equation \eqref{Pblm strong limit Darcy Omega 2} was already
shown in \eqref{Eq lower dimensional Darcy-type}, however the statements
\eqref{Pblm strong limit divergence Omega 2} and
\eqref{Pblm strong limit flux boundary condition} need further discussion.

\subsection{The interface integrals}\label{Sec interface integrals setting}

\begin{definition}\label{Def manifold spaces} \rm

Let $G$, $\zeta$ and $\Upsilon$ be as in definition
\ref{Def eligible surface and generated fissure}, define the spaces
\begin{subequations}\label{Def lower dimensional manifold isomorphisms}
\begin{gather}\label{Def L2 manifold}
L^{2}(\Upsilon) := \{h:\Upsilon\to\mathbb{R}: \int_{\Upsilon} h^{2}(s)  \,dS <+\infty  \} ,\\
\label{Def H1 manifold}
H^{1}(\Upsilon) := \{h  \in L^{2}(\Upsilon): \boldsymbol{\widetilde{\nabla}} h \in L^{2}(\Upsilon)
 \times L^{2}(\Upsilon)\} , \\
\label{Def H zero one manifold}
H^{1}_0(\Upsilon) := \{h\in H^{1}(\Upsilon): h|_{\partial G} = 0\} .
\end{gather}
\end{subequations}
Here $\boldsymbol{\widetilde{\nabla}}$ indicates the gradient with respect to the variables
$(x_{1}, x_{2})$ contained in $G$.
\end{definition}

We have the following isomorphism result.

\begin{theorem}\label{Th manifold isomorphism}
Let $\mathcal{O}$ and $\Upsilon$ be as in definition
\ref{Def eligible surface and generated fissure}.
Consider the natural embedding $\jmath: \Upsilon \to \mathcal{O}$
defined by $\jmath (\widetilde{\mathbf{x}}, \zeta(\widetilde{\mathbf{x}})) := \widetilde{\mathbf{x}}$ and the map
\begin{equation}\label{Def embedding}
\varphi \mapsto \varphi \circ \jmath
\end{equation}
Then
\begin{itemize}
\item[(i)]
The embedding \eqref{Def embedding} is an isomorphism between
$L^{2}(\mathcal{O})$ and $L^{2}(\Upsilon)$;

\item[(ii)]
The embedding \eqref{Def embedding} is an isomorphism between
$H^{1}(\mathcal{O})$ and $H^{1}(\Upsilon)$;

\item[(iii)]
The embedding \eqref{Def embedding} is an isomorphism between
$H_0^{1}(\mathcal{O})$ and $H^{1}_0(\Upsilon)$.
\end{itemize}
\end{theorem}

\begin{proof}
By definition $\jmath:  \Upsilon \to \mathcal{O}$ is linear and bijective,
therefore the map \eqref{Def embedding} is bijective between spaces of functions.

(i) By the hypothesis that $\zeta$ satisfies
$C_{1} = \operatorname{ess\,inf}\{\boldsymbol{\widehat{n}}(s)\cdot\boldsymbol{\widehat{k}}: s\in \Upsilon\} > 0$ then,
for any $\phi \in L^{2}(\Upsilon)$
\begin{equation*}%\label{Def embedding}
\int_{\Upsilon}  (\phi\circ \jmath)^{2} \,dS =
\int_{\mathcal{O}} (\boldsymbol{\widehat{n}}(\widetilde{\mathbf{x}})\cdot\boldsymbol{\widehat{k}})^{-1}  \phi^{2}(\widetilde{\mathbf{x}}) d\widetilde{\mathbf{x}}\leq
\frac{1}{C_{1}}\int_{\mathcal{O}}  \phi^{2}(\widetilde{\mathbf{x}}) d\widetilde{\mathbf{x}} .
\end{equation*}
The inequality above gives the continuity of the application
$\varphi\mapsto \varphi\circ j$. By Banach's inversion theorem the map
is an isomorphism.

(ii) By definition $\boldsymbol{\widetilde{\nabla}} (\phi \circ \jmath)  = \boldsymbol{\widetilde{\nabla}} \phi (\widetilde{\mathbf{x}})$
 holds for any $\phi\in H^{1}(\mathcal{O})$.

(iii) Is immediate from (ii).
\end{proof}


Choose $q\in Q_0$ supported inside $\Omega_{2}$ and test  equation
\eqref{Pblm lower dimensional variational statement 2}; hence
\begin{equation} \label{Eq tested divergence}
-\int_{\Omega_{2}}  \mathbf{\widetilde{u}}^2\cdot\boldsymbol{\widetilde{\nabla}} q
- \int_{\Gamma_1}  (\mathbf{u}^1\cdot\boldsymbol{\widehat{n}} )  q^2\,dS
+ \int_{\Gamma_{2}}  (\mathbf{u}^1\cdot\boldsymbol{\widehat{n}} )  q^2\,dS
= \int_{\Gamma_1\cup \Gamma_{2}}  f_{\Gamma_1}  q^2 \,dS .
\end{equation}
We focus on the first term of the left-hand side. First $\partial_{z} q^2 = 0$
implies $\mathbf{\widetilde{u}}^2\cdot\boldsymbol{\widetilde{\nabla}}  q^2 = \mathbf{u}^2\cdot\boldsymbol{\nabla}  q^2$, then
\begin{equation}\label{Eq integration by parts}
- \int_{\Omega_{2}}\mathbf{u}^2\cdot\boldsymbol{\nabla} q
= \int_{\Omega_{2}} \boldsymbol{\nabla}\cdot \mathbf{u}^2   q^2  -
\int_{\partial\Omega_{2}}  q^2 \mathbf{u}^2\cdot \boldsymbol{\widehat{\nu}}\, dS .
\end{equation}
The two summands on the right-hand side are treated separately.
For the first summand the independence from the variable $z$ implies
that $\boldsymbol{\nabla}\cdot \mathbf{u}^2 = \boldsymbol{\widetilde{\nabla}} \cdot \mathbf{\widetilde{u}}^2$. The fact $\mathbf{u}^2\cdot\boldsymbol{\widehat{n}} = 0$
in $\Omega_{2}$ gives $\mathbf{\widetilde{u}}^2 =  M^{ T,\boldsymbol{\tau}} \mathbf{u}_{\tau}^2$. Thus
\begin{equation}\label{Eq reduction of divergence term}
\int_{\Omega_{2}}  \boldsymbol{\widetilde{\nabla}} \cdot \mathbf{\widetilde{u}}^2   q^2\,d\mathbf{x} %\\
= h \int_{G}  \boldsymbol{\widetilde{\nabla}}\cdot ( M^{ T,\boldsymbol{\tau}}\mathbf{u}_{\tau}^2)   q^2\,d\widetilde{\mathbf{x}}
= h \int_{\Gamma_1}   (\boldsymbol{\widehat{n}}\cdot\boldsymbol{\widehat{k}}) \boldsymbol{\widetilde{\nabla}}\cdot ( M^{ T,\boldsymbol{\tau}}\mathbf{u}_{\tau}^2)   q^2\,dS .
\end{equation}
The boundary term in \eqref{Eq integration by parts} can be written as
\begin{equation*}%\label{Def embedding}
-\int_{\Gamma_1\cup \Gamma_{2}}  q^2 \mathbf{u}^2\cdot \boldsymbol{\widehat{\nu}}\,dS
-\int_{\partial\Omega_{2} - (\Gamma_1 \cup\Gamma_{2}) }  q^2 \mathbf{u}^2\cdot \boldsymbol{\widehat{\nu}}\,dS .
\end{equation*}
The first summand vanishes since $\mathbf{u}^2\cdot \boldsymbol{\widehat{n}} = 0$ in $\Omega_{2}$.
The boundary piece described in the second summand is a vertical wall,
then $\boldsymbol{\widehat{\nu}}\cdot\boldsymbol{\widehat{k}} = 0$ and it can be identified with the outwards
normal vector to the set $G\subseteq\mathbb{R}^{2}$. Moreover, from the
independence of the integrand with respect to the variable $z$, the surface
integral can be collapsed to a line integral over $\partial G$.
Combining these observations with \eqref{Eq reduction of divergence term}
and \eqref{Eq integration by parts}, the equation \eqref{Eq tested divergence}
transforms into
\begin{align*}
& h \int_{\Gamma_1} (\boldsymbol{\widehat{n}}\cdot\boldsymbol{\widehat{k}}) \boldsymbol{\widetilde{\nabla}}\cdot ( M^{ T,\boldsymbol{\tau}}\mathbf{u}_{\tau}^2)   q^2\,dS
-  h \int_{\partial G}  q^2 \mathbf{u}\cdot \boldsymbol{\widehat{\nu}}\,dC\\
&- \int_{\Gamma_1}(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}} )  q^2\,dS
+ \int_{\Gamma_{2}}(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}} )  q^2\,dS\\
&= \int_{\Gamma_1\cup \Gamma_{2}}f_{\Gamma}  q^2 dS,
\end{align*}
where $d C$ is the arc-length measure on $\partial G$.
The isomorphisms provided by theorem \ref{Th manifold isomorphism}
imply that the quantifier $ q^2|_{\Gamma_1} $ can hit any function in
the space $ H_0^{1}(\Gamma_1)$.  Therefore, the equation
\eqref{Pblm strong limit divergence Omega 2} follows. Finally, using again
theorem \ref{Th manifold isomorphism}, the trace of the test function
$ q^2|_{\Gamma_1} $ can hit any function in the space $ H_0^{1}(\Gamma_1)$,
and combined with equation \eqref{Pblm strong limit divergence Omega 2}
we obtain \eqref{Pblm strong limit flux boundary condition}.

\subsection{Strong convergence of solutions}\label{Sec strong convergence}

\begin{theorem}\label{Th strong convergence of solutions}
Under the hypotheses
\begin{equation} \label{Hyp strong convergence of forcing terms}
 \| F^{ \epsilon, 1}-F^{ 1}\| _{0, \Omega_{1}}\to 0 , \quad
\| f^{ \epsilon}_{\Gamma} - f_{\Gamma}\| _{0, \Gamma}\to 0, \quad
\| \mathbf{g}^{ \epsilon}- \mathbf{g}\| _{0, \Omega}\to 0 ,
\end{equation}
the solutions $\mathbf{u}^{\epsilon}, p^{\epsilon}$ satisfy the following strong convergence
statements
\begin{equation}\label{Hyp strong convergence of solutions}
\begin{gathered}
\| \mathbf{u}^{\epsilon,1} - \mathbf{u}^1\| _{0, \Omega_{1}}\to 0   , \quad
\| \epsilon \mathbf{u}^{\epsilon,2} - \mathbf{u}^2 \|_{0 , \Omega_{2}}\to 0   , \\
 \| p^{\epsilon,1} -  p^1\|_{1, \Omega_{1}}\to 0   , \quad
 \| p^{\epsilon,2} -  p^2\|_{1, \Omega_{2}}\to 0 .
\end{gathered}
\end{equation}
\end{theorem}

The proof of the above theorem uses exactly the same arguments presented in
\cite[Theorem 3.2]{MoralesShow2}, and it is omitted.

Finally, assume that $\mathbf{u}^2_{\tau}\neq 0$ and consider the quotients:
\begin{equation} \label{Eq blow up relation}
\frac{\| \mathbf{u}^{\epsilon,2}_{\tau}  \|_{0,\Omega_2}}{\| \mathbf{u}^{\epsilon,2}\cdot\boldsymbol{\widehat{n}} \|_{0, \Omega_2}}=
\frac{\| \epsilon \mathbf{u}^{\epsilon,2}_{\tau} \|_{0, \Omega_2}}
{\| \epsilon \mathbf{u}^{\epsilon,2}\cdot\boldsymbol{\widehat{n}} \|_{0, \Omega_2}}>
\frac{\| \mathbf{u}^2_{\tau} \|_{0, \Omega_2}-\delta}
{\| \epsilon \mathbf{u}^{\epsilon,2}\cdot\boldsymbol{\widehat{n}} \|_{0, \Omega_2}}>0\,.
\end{equation}
The lower bound holds true for $\epsilon>0$ small enough and adequate $\delta>0$.
Therefore, we conclude that the magnitudes' ratio of the tangential over
the normal components of the flux blows-up to infinity; i.e., the flow
in the thin channel is predominantly tangential. Finally if $\mathbf{u}^2_{\tau}=0$,
unlike the analysis for flat interfaces presented in \cite{MoralesShow2},
no conclusions can be obtained because of the complexity introduced by the geometry
of the fissure.

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.45\textwidth]{fig3}
\end{center}  % LowerDimensionalSingleFissure.pdf
\caption{System of 2-D manifold fissures}
\label{Fig Lower dimensional Fissured Medium}
\end{figure}

\section{A problem with two dimensional manifolds}
\label{Sec lower dimensional mixed formulation}

In this section, using the independence of the limit functions with respect
to $z$ in $\Omega_{2}$, it will be shown that the limiting problem
\eqref{Pblm limit problem strong form} can be formulated as a system
coupling Darcy flow in three dimensions, with tangential flow hosted in a
 piecewise $C^{1}$ surface, as depicted in figure
\ref{Fig Lower dimensional Fissured Medium}. First we introduce the geometry,
recall that
$\Theta_{2} - h  \boldsymbol{\widehat{k}} = \{\omega - k  \boldsymbol{\widehat{k}}: \omega \in \Omega_{2} \}$
and consider the domain
\begin{equation}\label{Def fissured medium region l-d}
\vartheta:=\Theta_{1}\cup (\Theta_{2} - h  \boldsymbol{\widehat{k}}),\quad
\vartheta_{FR}  :=  \vartheta\cup\Gamma_{1} .
\end{equation}
Additionally we introduce the notation $\Gamma_{1}^{+}$,
 $\Gamma_{1}^{-}$ for the upper and lower faces of the piecewise surface $\Gamma_{1}$.

\subsection{Spaces of functions and isomorphisms}
\label{Sec l-d spaces of functions and morphisms}

\begin{definition}\label{Def l-d fissured system spaces} \rm
We define the following spaces for velocity and pressure
\begin{subequations}\label{Eq l-d fissured system spaces}
\begin{gather}
\begin{aligned}\label{Def l-d fissured space of velocities}
\mathbf{V}_{f} := \big\{&\mathbf{v}\in \mathbf{L}^2(\vartheta_{FR}):\boldsymbol{\nabla}\cdot
\mathbf{v}^1 |_{\Theta_{1}}\in \mathbf{L}^2(\Theta_{1}), \boldsymbol{\nabla}\cdot \mathbf{v}^1 |_{\Theta_{2} - h
\boldsymbol{\widehat{k}}}\in \mathbf{L}^2(\Theta_{2} - h \boldsymbol{\widehat{k}}) ,\\
&\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{+}} ,  \mathbf{v}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{-}} \in L^{2}(\Gamma_1) ,
\mathbf{v}^2|_{\Gamma_1}\in \mathbf{L}^2(\Gamma_1)
\big\} ,
\end{aligned}\\
\label{Def l-d fissured space of pressures}
Q_{f} := \{q\in L^{2}(\vartheta_{FR}):  q|_{\Gamma_1}\in H^{1}(\Gamma_1)\} .
\end{gather}
Endowed with the norms coming from the natural inner products
\begin{gather} \label{Def l-d norm fissured space of velocities}
\begin{aligned}
\|\mathbf{v}\|_{\mathbf{V}_{f}} &:= \big\{\| \mathbf{v} \|_{\mathbf{L}^2(\vartheta_{FR})}^{ 2}
+\| \boldsymbol{\nabla}\cdot \mathbf{v} \|_{L^{2}(\vartheta_{FR})}^{ 2}
+\| \mathbf{v}\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{+}} \|_{L^{2}(\Gamma_1)}^{ 2}\\
&\quad +
\| \mathbf{v}\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{-}} \|_{L^{2}(\Gamma_1)}^{ 2}+
\| \mathbf{v} \|_{\mathbf{L}^2(\Gamma_1)}^{ 2}\big\}^{1/2} ,
\end{aligned}\\
\label{Def l-d norm fissured space of pressures}
\| q \|_{Q_{f}} := \{\|  q \|_{L^{2}(\vartheta_{FR})}^{ 2}
+ \| q \|_{H^{1}(\Gamma_1)}^{ 2}\}^{1/2} .
\end{gather}
\end{subequations}
\end{definition}

\begin{remark}\label{Rem structure of the collapsed space} \rm
Note that definition \eqref{Def l-d fissured space of velocities} requires
only $\mathbf{v}\in H_{\rm div}(\Theta_{1})$ and $\mathbf{v}\in H_{\rm div}(\Theta_{2}- h \boldsymbol{\widehat{k}})$; i.e.,
 the divergence is a square integrable function only on these subdomains.
Therefore, both normal traces $\mathbf{v}\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{+}}$ and
$\mathbf{v}\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{-}}$ make sense in $H^{-1/2}(\Gamma_1)$, but we require
the extra condition of been in $L^{2}(\Gamma_1)$. We do not demand the global
 condition $\mathbf{v}\in H_{\rm div}(\vartheta_{FR})$ because this would imply the continuity
of the normal traces across a surface; i.e.,
$\mathbf{u}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{+}} =  \mathbf{u}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{-}}$. Such condition
can not model jumps across the fissures as the normal stress balance
interface \eqref{Pblm strong normal stress balance} and the limit
 equation \eqref{Pblm strong limit divergence Omega 2}.
\end{remark}

Next define a change of variable based on piecewise rigid translations

\begin{definition}\label{Def translation change of variable} \rm
Let $\mathbf{x} = (\widetilde{\mathbf{x}}, x_{3}) $ and define the map $T:\Omega\to\mathbb{R}^{ 3}$
\begin{equation}\label{Eq translation change of variable}
T \mathbf{x}:= (\widetilde{\mathbf{x}}, x_{3} )\boldsymbol{\mathbbm{1}}_{\Theta_{1}}(\widetilde{\mathbf{x}},  x_3)
+ (\widetilde{\mathbf{x}}, x_{3} -  h)\boldsymbol{\mathbbm{1}}_{\Theta_{2}}(\widetilde{\mathbf{x}},  x_3)
+
(\widetilde{\mathbf{x}}, \zeta(\widetilde{\mathbf{x}}) ) \boldsymbol{\mathbbm{1}}_{\Theta( h, \Gamma_1)}(\widetilde{\mathbf{x}},  x_3)
\end{equation}
\end{definition}

\begin{theorem}\label{Th l-d functions space isomorphisms}
\begin{itemize}
\item[(i)] The application $\mathbf{v} \mapsto \mathbf{v}\circ T$ is an isometric isomorphism from
$\mathbf{V}_0$ to $\mathbf{V}_{f}$.

\item[(ii)] The application $q\mapsto q\circ T$ is an isometric isomorphism
from $Q_0$ to $Q_{f}$.
\end{itemize}
\end{theorem}

\begin{proof}
(i) The proof is a direct application of part (i) in theorem
\ref{Th manifold isomorphism}. The only detail that needs further
clarification is to observe that
 \begin{gather*}
\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1}\mapsto (\mathbf{v}^1\circ T)\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{-}}, \\
\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_{2}} =
\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1+ h}\mapsto (\mathbf{v}^1\circ T)\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{+}}\,.
\end{gather*}

 (ii) It is a direct application of parts (i) and (ii) in
theorem \ref{Th manifold isomorphism}.
\end{proof}

\subsection{The lower dimensional mixed problem}\label{Sec l-d mixed problem}

By the previous theorem, problem \eqref{Pblm limit system}
is equivalent to the following mixed problem with a piecewise $C^{1}$
coupling interface: Find $p\in Q_{f}$ and $\mathbf{u}\in \mathbf{V}_{f}$ such that
\begin{subequations}\label{Pblm reduced dimension variational statement}
\begin{gather}\label{Pblm reduced dimension variational statement 1}
\begin{aligned}
&\int_{\vartheta} a_1  \mathbf{u} \cdot \mathbf{v} - \int_{\vartheta} p  \boldsymbol{\nabla}\cdot\mathbf{v}
+  h\int_{\Gamma_1} (\boldsymbol{\widehat{n}}\cdot\boldsymbol{\widehat{k}} )(a_2  \mathbf{u}_{\tau}^2 + ( M^{ T,\boldsymbol{\tau}})'
\widetilde{\boldsymbol{\nabla}} p +\mathbf{g}_{\tau})\cdot \mathbf{v}_{\tau}^2  \,dS\\
&+\alpha\int_{\Gamma_1}\big[(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{+}})(\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{+}})
+(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{-}})(\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{-}})\big] \,dS\\
&- \int_{\Gamma_1} p^2 \big[(\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{+}})
- (\mathbf{v}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{-}})\big] \,dS\\
&= - \int_{\vartheta}\mathbf{g} \cdot \mathbf{v} - \int_{\Gamma_{1}}\mathbf{g}_{\tau} \cdot \mathbf{v}_{\tau}^2 \,dS\,,
\end{aligned}\\
\label{Pblm reduced dimension variational statement 2}
\begin{aligned}
&\int_{\vartheta}\boldsymbol{\nabla}\cdot\mathbf{u} q
-    h \int_{\Gamma_1} (\boldsymbol{\widehat{n}}\cdot\boldsymbol{\widehat{k}} )
 M^{ T,\boldsymbol{\tau}}\mathbf{u}_{\tau}^2 \cdot \boldsymbol{\widetilde{\nabla}}  q^2  \,dS
+ \int_{\Gamma_1}\big[(\mathbf{u}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{+}})
- (\mathbf{u}^1\cdot\boldsymbol{\widehat{n}}|_{\Gamma_1^{-}})\big]  q^2 \,dS \\
&= \int_{\vartheta}F^{  1} q
+ \int_{\Gamma_{1}} f_{\Gamma}   q^2 \,dS\quad
\text{for all } q\in Q_{f} ,\; \mathbf{v}\in \mathbf{V}_{f} .
\end{aligned}
\end{gather}
\end{subequations}

Finally the equivalence of problems \eqref{Pblm limit system}
and \eqref{Pblm lower dimensional variational statement}  gives
the well-posedness of the system above.



\section{Final remarks and future work}
\label{Sec limitations and future work}

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.9\textwidth]{fig4}
\end{center}  % BothFissures.pdf
\caption{Translation generated fissures}\label{Fig Fissure with high gradients}
\end{figure}

(i) Giving the adequate definitions, the results can be generalized
immediately to a system of multiple fissures, such as the one depicted
in figure \ref{Fig Fractured Medimu}. The formulation presented in this work
can handle large amounts of information in a remarkably efficient way.
One of the main reasons is the notation introduced by Showalter in
\cite{MoralesShow2} for the description of function spaces.

(ii) Our results can be generalized immediately to the $\mathbb{R}^{ N}$-setting
using the same arguments presented here. The problems have analogous structure.

(iii)
The approach based on analytic semigroups theory presented in section
\cite{MoralesShow2}, can be directly applied to this case, in order to
model the time dependent problem for totally fissured systems with singularities.

(iv) Although the mathematical analysis is solid, the approach used throughout
the paper stops been suitable for surfaces with high gradients, such as the one
depicted in the right-hand side of figure \ref{Fig Fissure with high gradients}
where $\boldsymbol{\widehat{n}}_{2}\cdot\boldsymbol{\widehat{k}} \ll \boldsymbol{\widehat{n}}_{1}\cdot\boldsymbol{\widehat{k}}$.
In this case, the translation in the direction $\boldsymbol{\widehat{k}}$ generates a fissure
whose cross section areas can be very different from one piece to another;
 i.e., $A_{2} \ll A_{1}$. Such a fissure is not realistic. On the other hand,
consider a fissure such as the one depicted in the left-hand side of
figure \ref{Fig Fissure with high gradients}. Here the translation is made
in the bisector vector direction
\[
\boldsymbol{\widehat{e}} \equiv \frac{1}{| \frac{\boldsymbol{\widehat{n}}_{1}+\boldsymbol{\widehat{n}}_{2}}{2} |} \frac{\boldsymbol{\widehat{n}}_{1}+\boldsymbol{\widehat{n}}_{2}}{2}
\]
This process generates a more realistic fissure.

(v) Demanding the fissures to be defined by the parallel translation of a
surface in a fixed direction, is definitely a step forward with respect to
previous achievements, however it is still a restrictive hypothesis for
modeling the phenomenon in natural geological formations.

(vi) Setting the problem in the mixed variational formulation used here,
can be easily extended to systems with fissures described by a very general
type of geometry. However, the difficulty of the asymptotic analysis,
for upscaling purposes, increases substantially.

(vii) Such questions will be addressed in future work by the introduction
of correction factors, obtained comparing the flow energy dissipation
in a real fissure and an artificial one e.g. replacing the presence of
the fissure in the left-hand side of Figure \ref{Fig Fissure with high gradients}
with the one on the right-hand side affected by a correction factor.
In the same way, fissures defined by walls which are not rigid translations
of a common surface, will be compared to a fissure generated by the vertical
translation of an ``average surface'' having the same ``average width''.


\subsection*{Acknowledgements}\label{Sec Acknowledgements}

The author thanks to Universidad Nacional de Colombia, Sede Medell\'in
for partially supporting this work under the projects HERMES 17194 and HERMES 14917
as well as the Department of Energy, Office of Science, USA for partially
supporting this work under grant 98089. Finally, the author wishes to thank
Professor Ralph E. Showalter from the Mathematics Department at Oregon
State University for his helpful insight, observations and sound suggestions.

\begin{thebibliography}{00}

\bibitem{Allaire2009}
Gr\'egorie Allaire, Marc Briane, Robert Brizzi, and Yves Capdeboscq;
\newblock Two asymptotic models for arrays of underground waste containers.
\newblock {\em Applied Analysis}, 88 (no. 10-11):1445--1467, 2009.

\bibitem{ArbBrunson2007}
Todd Arbogast and Dana Brunson.
\newblock A computational method for approximating a {D}arcy-{S}tokes system
  governing a vuggy porous medium.
\newblock {\em Computational Geosciences}, 11, No 3:207--218, 2007.

\bibitem{ArbLehr2006}
Todd Arbogast and Heather Lehr.
\newblock Homogenization of a darcy-stokes system modeling vuggy porous media.
\newblock {\em Computational Geosciences}, 10, No 3:291--302, 2006.

\bibitem{Gunzburger2009}
Nan Chen, Max Gunzburger, and Xiaoming Wang.
\newblock Asymptotic analysis of the differences between the {S}tokes-{D}arcy
  system with different interface conditions and the {S}tokes-{B}rinkman
  system.
\newblock {\em Journal of Mathematical Analysis and Applications}, 368
  (2):658--676, 2009.

\bibitem{Gatica2009}
Gabriel~N. Gatica, Salim Meddahi, and Ricardo Oyarz\'ua.
\newblock A conforming mixed finite-element method for the coupling of fluid
  flow with porous media flow.
\newblock {\em IMA Journal of Numerical Analysis}, 29, 1:86--108, 2009.

\bibitem{GiraultRaviartFEM}
V.~Girault and P.-A. Raviart.
\newblock {\em Finite element approximation of the {N}avier-{S}tokes
  equations}, volume 749 of {\em Lecture Notes in Mathematics}.
\newblock Springer-Verlag, Berlin, 1979.

\bibitem{Hornung}
Ulrich Hornung, editor.
\newblock {\em Homogenization and Porous Media, Ulrich Hornung editor},
  volume~6 of {\em Interdisciplinary Applied Mathematics}.
\newblock Springer-Verlag, New York, 1997.

\bibitem{ZhaoQin}
ZhaoQin Huang, Jun Yao, YaJun Li, ChenChen Wang, and XinRui L\"{u}.
\newblock Permeability analysis of fractured vuggy porous media based on
  homogenization theory.
\newblock {\em Science China Technological Sciences}, 53 (3):839--847, 2010.

\bibitem{Yotov}
W.~J. Layton, F.~Scheiweck, and I.~Yotov.
\newblock Coupling fluid flow with porous media flow.
\newblock {\em SIAM Journal of Numerical Analysis}, 40 (6):2195--2218, 2003.

\bibitem{Levy83}
Th{\'e}r{\`e}se L{\'e}vy.
\newblock Fluid flow through an array of fixed particles.
\newblock {\em International Journal of Engineering Science}, 21:11--23, 1983.

\bibitem{Loredana1}
J.~San Mart{\'i}n, J.-F. Scheid, and L.~Smaranda.
\newblock A modified lagrange-galerkin method for a fluid-rigid system with
  discontinuous density.
\newblock {\em Numerische Mathematik}, 122 (2):341--382, 2012.

\bibitem{JaffRob05}
Vincent Martin, J{\'e}r{\^o}me Jaffr{\'e}, and Jean~E. Roberts.
\newblock Modeling fractures and barriers as interfaces for flow in porous
  media.
\newblock {\em SIAM J. Sci. Comput.}, 26(5):1667--1691, 2005.

\bibitem{Mikelic89}
A.~Mikeli{\'c}.
\newblock A convergence theorem for homogenization of two-phase miscible flow
  through fractured reservoirs with uniform fracture distribution.
\newblock {\em Applicable Anal.}, 33:203--214, 1089.

\bibitem{Morales1}
Fernando A. Morales.
\newblock Analysis of a coupled Darcy multiple scale flow model under geometric perturbations of the interface.
\newblock {\em Journal of Mathematics Research},
  5(4):11--25, 2013. DOI: 10.5539/jmr.v5n4p11

\bibitem{MoralesShow1}
Fernando Morales and Ralph Showalter.
\newblock The narrow fracture approximation by channeled flow.
\newblock {\em Journal of Mathematical Analysis and Applications},
  365:320--331, 2010.

\bibitem{MoralesShow2}
Fernando Morales and Ralph Showalter.
\newblock Interface approximation of {D}arcy flow in a narrow channel.
\newblock {\em Mathematical Methods in the Applied Sciences}, 35:182--195,
  2012.

\bibitem{SP80}
Enrique S\'anchez-Palencia.
\newblock {\em Nonhomogeneous media and vibration theory}, volume 127 of {\em
  Lecture Notes in Physics}.
\newblock Springer-Verlag, Berlin, 1980.

\bibitem{Showalter77}
R.~E. Showalter.
\newblock {\em Hilbert space methods for partial differential equations},
  volume~1 of {\em Monographs and Studies in Mathematics}.
\newblock Pitman, London-San Francisco, CA-Melbourne, 1977.

\bibitem{Showalter97}
R.~E. Showalter.
\newblock {\em Monotone operators in {B}anach space and nonlinear partial
  differential equations}, volume~49 of {\em Mathematical Surveys and
  Monographs}.
\newblock American Mathematical Society, Providence, RI, 1997.

\bibitem{Show97}
R.E. Showalter.
\newblock {\em Microstructure Models of Porous Media. In Ulrich Hornung editor
  Homogenization and Porous Media}, volume~6 of {\em Interdisciplinary Applied
  Mathematics}.
\newblock Springer-Verlag, New York, 1997.

\end{thebibliography}


\end{document}