\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2007(2007), No. 147, pp. 1--30.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2007 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2007/147\hfil Multiscale systems] {Multiscale elliptic-parabolic systems \\ for flow and transport} \author[M. Peszy{\'n}ska, R. E. Showalter \hfil EJDE-2007/147\hfilneg] {Ma{\l}gorzata Peszy{\'n}ska, Ralph E. Showalter} \address{Ma{\l}gorzata Peszy{\'n}ska \newline Department of Mathematics \\ Oregon State University, Corvallis, OR 97331, USA} \email{mpesz@math.oregonstate.edu} \address{Ralph E. Showalter \newline Department of Mathematics \\ Oregon State University, Corvallis, OR 97331, USA} \email{show@math.oregonstate.edu} \thanks{Submitted October 5, 2007. Published November 5, 2007.} \thanks{The first author was supported by grants 0511190 from the National Science Foundation, \hfill\break\indent and 98089 from Department of Energy, Office of Science.} \thanks{The second author was supported by grants 98089 and 9001997 from Department of Energy, \hfill\break\indent Office of Science.} \subjclass[2000]{76S05, 35B27, 74Q15, 35R10} \keywords {Upscaled porous media; double porosity models; \hfill\break\indent multiscale flow and transport; nonlocal dispersion} \begin{abstract} An upscaled elliptic-parabolic system of partial differential equations describing the multiscale flow of a single-phase incompressible fluid and transport of a dissolved chemical by advection and diffusion through a heterogeneous porous medium is developed without the usual assumptions of scale separation. After a review of homogenization results for the traditional low contrast and the $\epsilon^2$-scaled high contrast cases, the new discrete upscaled model based on local affine approximations is constructed. The resulting model is mass conserving and contains the effects of local advective transport as well as diffusion, it includes non-Fickian models of dispersion and works over a broad range of contrast cases. \end{abstract} \maketitle \allowdisplaybreaks \numberwithin{equation}{section} \newtheorem{lemma}{Lemma}[section] \newtheorem{corollary}[lemma]{Corollary} \newtheorem{proposition}[lemma]{Proposition} \newtheorem{remark}[lemma]{Remark} \newcommand{\dtt}[1]{\frac{\partial #1}{\partial t}} \newcommand{\ddt}[1]{\frac{d #1}{d t}} \newcommand{\lave}[2]{\langle #1 \rangle_{#2}} \newcommand{\ivec}[1]{\left(#1\right)_{i=1}^{N_{\rm incl}}} \newcommand{\abs}[1]{|#1|} \newcommand{\norm}[2]{\|#1\|_{#2}} \newcommand{\snorm}[2]{|#1|_{#2}} \newcommand{\mpcite}[2]{{[\cite{#1}, {\it #2}]}} \section{Introduction} \label{sec-prelim} We construct a new family of differential models of flow and transport in heterogeneous porous media, a system consisting of an elliptic equation for stationary flow coupled to a parabolic equation for the transient advection-diffusion. The new features of these models are (i) they work over a range of coefficient and geometry scales, (ii) their discrete form is amenable to numerical simulation, (iii) the mass conservation property is preserved in the process of upscaling, and (iv) they retain the variational structure and well-posedness of the initial-boundary-value problem. Our work was originally motivated by a particular experiment \cite{Hagg04} which provides a unique testbed for multiscale modeling in the presence of not well separated scales. However, the results immediately apply to a general class of physical phenomena of flow and transport in variously heterogeneous geological formations. The separation of scales is not assumed in this paper. By this we understand any of the following. First, the ratio $\epsilon$ of the diameter of a typical cell or representative elementary volume to the diameter of the physical domain is very small and in traditional homogenization or volume averaging techniques is driven to $0$ to obtain the upscaled model. By contrast, here we consider a fixed $\epsilon=\epsilon_0>0$. Second, the separation of scales of flow or diffusion requires the ratio of fast to slow coefficients to be either bounded, the low contrast case, or $\mathcal{O}(\epsilon^{-2})$, the high contrast case. The latter leads to the double porosity models. However, here we permit a wide range of such ratios from low through intermediate to high constrast of coefficients. Thirdly, with strongly separated scales the dominating effect is diffusion, whereas here we capture a range of advection--diffusion--dispersion phenomena. In addition to being computationally tractable, these new model systems give qualitative information about macroscopically observable quantities and their time scales and rates, and they allow us to observe the transition between various regimes of flow and transport. Traditional methods of homogenization \cite{BLP,SP,JKO,Allaire92} strive to determine coefficients of an effective partial differential equation (PDE). When the fine-scale solutions vary substantially both in space and time, the solutions of these effective equations fail to convey essential information. The double-porosity models \cite{BZK60,Vogt,Arb88,Arb89,Arb89r,ADH,DougArb,HornSho90,ShowWalk91a} retain nonlocal effects in time from cell problems on the original scale. In hydrology applications related models, also called {\em multi-rate} models or {\em power-law} models, are widely used; see \cite{Hagg95,Hagg00,roy-tailing,roy-hyporheic}. The goal is to capture both short-term and long-term tailing due to diffusive storage and adsorption. The need for new models has been pointed out for example in \cite{Russ88,ClarkShow94}, in \mpcite{Arbogast97}{p.217} and in \cite{BP98} where the regimes were studied carefully. The recently proposed models in \cite{Panfilov} and \cite{Hagg04,wood.cmg} work well for some regimes of flow and transport but not across all. We describe the model problem below, a heterogeneous system with combined fast and slow flow regimes. In Section~\ref{sec-review} we review known homogenization results for the elliptic and parabolic subproblems of that coupled system. We shall exploit this theory but will not contribute to it. In Section~\ref{sec-discrete} we propose the differential system of central interest in this paper, one PDE coupled across interfaces to a collection of cell problems. Local affine approximations extend the traditional constant approximations of constraints on cell interfaces and introduce additional nonlocal terms which are carefully constructed so that the mass is preserved. The variational setting leads to well-posedness of the system. In Section~\ref{sec-convolution} we compute the nonlocal terms as convolutions and derive the continuous limit of the generic effective differential model. The model has a form strikingly similar to those discussed in \cite{GoltzRoberts87,Cushman93,ELW00,SanCarr04} and displays the various flow regimes \cite{BP98,Panfilov}. Finally, in Section~\ref{sec-newmodel}, we return to the original flow and transport system and construct the effective model. This includes the velocity and dispersion terms which are lost by traditional piecewise constant interface approximations. Throughout the paper we use the following notation. Let $D \subset \mathbb{R}^2$ be an open bounded set with boundary $\partial D$. (We restrict ourselves to $\mathbb{R}^2$ for simplicity only; this is consistent with \cite{Hagg04}). In area integrals we use the symbol $dA$ and use $dS$ in surface integrals. The characteristic function of a set $D$ is $\chi_D$, and $ \lave{f}{D} \equiv \frac{1}{\abs{D}} \left( \int_{D} f(\mathbf{x})\,dA \right)$ is the average of $f$ over $D$. At times we use the notion of a translate $D(\mathbf{x})$ with centroid at the point $\mathbf{x}$. In such cases the spatial variable is denoted by $\mathbf{y} \in D(\mathbf{x})$. \section*{The Model Problem} Consider flow and transport in a heterogeneous porous medium, an open bounded domain $\Omega \subset \mathbb{R}^2$. We denote by $\mathbf{n}$ the unit normal vector out of $\Omega$ on the boundary $\partial \Omega$. The flow of water is described by conservation of mass and Darcy's law, respectively, % \begin{equation} \label{eq-flow} \nabla \cdot \mathbf{v} = 0, \quad \mathbf{v} = - \mathbf{K} \nabla p, \; \mathbf{x} \in \Omega\,, \end{equation} % an elliptic equation for the pressure $p(\mathbf{x})$, where $\mathbf{v}=(v_1,v_2)=[v_1,v_2]^T$ is the volumetric flux (velocity), and $\mathbf{K}:\mathbb{R}^2 \mapsto \mathbb{R}^{2 \times 2}$ is the symmetric uniformly positive definite conductivity representing permeability divided by fluid viscosity. The flowing water contains a dissolved conservative dye of concentration $c(\mathbf{x},t)$. The associated model of transient advection--diffusion--dispersion is the parabolic equation \cite{Bear72,Ewing83} %% \begin{eqnarray} \phi \frac{\partial c}{\partial t}+ \nabla \cdot ( \mathbf{v} c - \mathbf{D}(\mathbf{v}) \nabla c) = 0, \quad \mathbf{x} \in \Omega\,, \label{eq-transport} \end{eqnarray} %% The coefficient $\phi$ is the porosity of the medium, and the diffusion-dispersion tensor is given by % \begin{eqnarray} \mathbf{D} = \mathbf{D}(\mathbf{v}) \equiv d_{\rm mol} \mathbf{I} + \snorm{\mathbf{v}}{} (d_{l} \mathbf{E}(\mathbf{v}) + d_t (\mathbf{I} - \mathbf{E}(\mathbf{v}))). \label{eq-dispersion} \end{eqnarray} % Here $d_{mol},d_l,d_t$ are coefficients of molecular diffusivity, longitudinal and transversal dispersivity, respectively, and the dispersion tensor $\mathbf{E}(\mathbf{v}) = \frac{1}{\snorm{\mathbf{v}}{}^2} v_i v_j$ is a rank two tensor. The coupled problem \eqref{eq-flow}--\eqref{eq-transport} with appropriate boundary and initial conditions is well understood under standard assumptions of mild variation on coefficients. The difficulties arise when the coefficients $\phi,\mathbf{K}$ and, consequently, $\mathbf{v}$ and $\mathbf{D}(\mathbf{v})$, have a highly heterogeneous character and differ by several orders of magnitude. As in \cite{Hagg04}, we assume the coefficients are locally piecewise constant on an associated collection of local interfaces separating two disjoint regimes of flow; the subscripts $f,s$ are associated with the {\em fast}, and {\em slow} regions, $\Omega_f,\Omega_s$, respectively. These are disjoint open sets covering $\Omega$, $\overline{\Omega} = \overline{\Omega_f} \cup \overline{\Omega_s}$, with an interface $\Gamma_{fs} \equiv \partial \Omega_f \cap \partial \Omega_s$. The region $\Omega_f$ is connected, but $\Omega_s = \bigcup_{i=1}^{N_{\rm incl}} \Omega_{is}$ where each {\em inclusion} $\Omega_{is}$ is a connected and simply connected region with $\Omega_{is} \cap \Omega_{js} =\emptyset, \; i \neq j$. We also denote the local interfaces by $\Gamma_i \equiv \partial \Omega_{is} \cap \partial \Omega_f$, so $\Gamma_{fs} = \bigcup_{i} \Gamma_i$. In what follows $\mathbf{n}_i$ denotes the unit normal to $\Gamma_i$ pointing out of $\Omega_{is}$. In the context of fractured (fissured) media \cite{DougArb,Arbogast97,HornSho90,P92} $\Omega_f$ is called a {\em fracture system} and $\Omega_s$ the {\em matrix} composed of {\em blocks} $\Omega_{is}$, a {\em totally fissured medium} \cite{DPS97}. We further assume that $\Omega$ can be covered by a union of rectangular subdomains $\Omega_i, i=1, \ldots N_{\rm incl}$ with each $\Omega_i$ containing exactly one inclusion $\Omega_{is}$. We denote by $\Omega_{if} = \Omega_i \cap \Omega_f$ the {fast} part surrounding $\Omega_{is}$ so that $\Omega_i = \Omega_{is} \cup \Omega_{if} \cup \Gamma_i$. Furthermore, we assume that each $\Omega_i$ is congruent to a generic cell $\Omega_0$ of size $\epsilon_0=diam(\Omega_0)$ and that each $\Omega_{is}$ is congruent to a generic $\Omega_{0s}$. It follows that the interfaces $\Gamma_i$ are congruent to $\Gamma_0$ and $\Omega_{if}$ is congruent to some $\Omega_{0f}$. Such assumptions of periodicity of $\Omega$ make the model amenable to homogenization but are not required for our development. Denote by $\mathbf{x}^C_i$ the centroid of $\Omega_{is}$ and by $\hat{\chi}_i(\mathbf{x})$ the characteristic function of the cell $\Omega_i$. Also, we denote by $\theta_f=\frac{\abs{\Omega_{0f}}}{\abs{\Omega_0}}$ the volume fraction occupied by the ``fast'' region. Analogously $\theta_s=\frac{\abs{\Omega_{0s}}}{\abs{\Omega_0}}$ is defined with $\theta_s=1-\theta_f$. We consider scalar coefficients of the form %% \begin{eqnarray} \label{eq-pc} a(\mathbf{x}) = a(a_f,a_s; \mathbf{x}) \equiv a_f \chi_{\Omega_f} (\mathbf{x}) + a_s \chi_{\Omega_s}(\mathbf{x}) =\left\{ \begin{array}{ll} a_f, \; \mathbf{x} \in \Omega_f \\ a_s,\; \mathbf{x} \in \Omega_s \end{array}\right. , \; \mathbf{x} \in \Omega \end{eqnarray} %% where $a_f,a_s >0$ are given constants, as well as isotropic tensor valued coefficients %% \begin{eqnarray} \label{eq-pcti} \mathbf{A}(\mathbf{x}) \equiv A(A_f,A_s; \mathbf{x}) \mathbf{I} = A_f \chi_{\Omega_f} (\mathbf{x})\mathbf{I} + A_s \chi_{\Omega_s}(\mathbf{x}) \mathbf{I}, \; \mathbf{x} \in \Omega. \end{eqnarray} %% More general forms of \eqref{eq-pc}, \eqref{eq-pcti} taking different values in each inclusion $\Omega_{is}$ are, respectively, %% \begin{gather} \label{eq-pci} a(\mathbf{x}) \equiv a(a_f,\left(a_i\right)_{i=1}^{N_{\rm incl}}; \mathbf{x}) = a_f \chi_{\Omega_f} (\mathbf{x}) + \sum_i a_i \chi_{\Omega_{is}}(\mathbf{x}), \; \mathbf{x} \in \Omega\,, \\ \label{eq-pctii} \mathbf{A}(\mathbf{x}) \equiv A(A_f,\ivec{A_i}; \mathbf{x}) \mathbf{I} = A_f \chi_{\Omega_f} (\mathbf{x})\mathbf{I} + \sum_i A_i \chi_{\Omega_{is}}(\mathbf{x}) \mathbf{I}, \; \mathbf{x} \in \Omega \,. \end{gather} %% In particular, as in the experimental setup in \cite{Hagg04}, we assume that the coefficients $\phi,\mathbf{K}$ satisfy \eqref{eq-pc}, and \eqref{eq-pcti}, respectively. On the other hand, we allow the diffusion--dispersion coefficient $\mathbf{D}(\mathbf{x})$ to depend, in general, on $i$ as in \eqref{eq-pctii}, because it depends on the local velocity $\mathbf{v}_i$ which in turn depends on $\mathbf{K}_i$ and on the local flow in $\Omega_{is}$. The {\em transmission form} of the coupled system \eqref{eq-flow}--\eqref{eq-transport} displays the heterogeneity of the coefficients as well as the conservation of momentum and mass across interfaces, %% \begin{subequations} \label{eq-adv-dif-flo-int} \begin{eqnarray} \nabla \cdot \mathbf{v}_f = 0,\ \mathbf{v}_f = - \mathbf{K}_f \nabla p_f, \quad x \in \Omega_f, \\ %\label{eq-flow-i} \nabla \cdot \mathbf{v}_i = 0,\ \mathbf{v}_i = - K_s \nabla p_i, \quad x \in \Omega_i, \;\; i=1,\ldots N_{\rm incl} \\ p_i = p_f, \;\; \mathbf{v}_i \cdot \mathbf{n} = \mathbf{v}_f \cdot \mathbf{n}, \;\; x \in \Gamma_i, \end{eqnarray} \end{subequations} % \begin{subequations} \label{eq-adv-dif-tra-int} \begin{eqnarray} \phi_{f} \frac{\partial c_{f}}{\partial t} - \nabla \cdot ( \mathbf{D}_{f} \nabla c_{f} - \mathbf{v}_f c_{f}) = 0, \;\; \mathbf{x} \in \Omega_{f}, %\label{eq-multi-f} \\ \phi_{i} \frac{\partial c_{i}}{\partial t} - \nabla \cdot ( \mathbf{D}_{i} \nabla c_{i} - \mathbf{v}_i c_{i}) = 0, \;\; \mathbf{x} \in \Omega_{i}, \;\; i=1,\ldots N_{\rm incl} %\label{} \\ c_i = c_f, \;\; (\mathbf{D}_f \nabla c_f - \mathbf{v}_f c_{f}) \cdot \mathbf{n} = (\mathbf{D}_i \nabla c_i - \mathbf{v}_i c_{i}) \cdot \mathbf{n}, \;\; \mathbf{x} \in \Gamma_i\,. \end{eqnarray} \end{subequations} % The system \eqref{eq-adv-dif-flo-int}, \eqref{eq-adv-dif-tra-int} is the {\em exact discrete model} for the problem from \cite{Hagg04}. Theoretically, it can be solved numerically, if enough computational resources are available. Also theoretically, the computed values of $p$ and $c$ could be verified pointwise thanks to currently available experimental and visualization techniques such as the imaging equipment used in \cite{Hagg04}. The discussion of \cite{Hagg04} centers around identification and fitting of coefficients $\mathbf{v},\ \mathbf{D}(\mathbf{v})$ to known upscaled versions of \eqref{eq-transport}. Depending on the ratio $\frac{K_f}{K_s}$, three distinct regimes of flow and transport are discussed. These are, respectively, $\frac{K_f}{K_s} = 6, 300, 1800$, and are called respectively, the {\em low contrast}, {\em intermediate}, and {\em high contrast} cases. It is concluded that this approach gives satisfactory results in the low and high contrast regimes where, respectively, \eqref{eq-transport} and its double-porosity modification are used, but that neither of these is satisfactory in the intermediate case. In the last case it is impossible to fit the observed breakthrough curves with available models. The differential model developed below works across all regimes of flow and transport from low to high constrast and includes the intermediate regime. \section{Homogenization of problems with discontinuous coefficients} \label{sec-review} We review traditional homogenization of elliptic and parabolic partial differential equations. Various approaches to upscaling of periodically varying discontinuous coefficients exist, and the two main classes can be distinguished by the use or not of $\epsilon^2$-scaling. Here we review these results and motivate the need to go beyond them. A novel and pivotal observation is that the ``obstacle'' problems \cite{JKO} provide a bridge between classical models without scaling \cite{BLP,SP,JKO} and double-porosity models with $\epsilon^2$-scaling \cite{DougArb,ADH}. The following applies to both the elliptic problem for the flow part \eqref{eq-adv-dif-flo-int} of the coupled system and the parabolic problem for the transport part \eqref{eq-adv-dif-tra-int}; the discussion follows \mpcite{SP}{II.5}, \mpcite{JKO}{Ch.3}. Consider the parabolic problem, % \begin{subequations} \label{eq-parabolic} \begin{eqnarray} \label{eq-parabolic-model} \phi\frac{\partial u}{\partial t}- \nabla \cdot (\mathbf{B} \nabla u) &=& f, \; \mathbf{x} \in \Omega, \\ \label{eq-parabolic-extbc} u &=& 0, \;\; \mathbf{x} \in \partial \Omega, \\ \label{eq-parabolic-initf} u(\mathbf{x},0) &=& u_{0}(\mathbf{x}), \;\; \mathbf{x} \in \Omega\,, \end{eqnarray} %% \end{subequations} % a special case of \eqref{eq-transport} in which dispersion and advection are ignored and with a source/sink term $f: \Omega \mapsto \mathbb{R}$. Assume that the coefficient $\phi (\mathbf{x})\equiv \phi(\phi_f,\phi_s;\mathbf{x})$ satisfies \eqref{eq-pc} and that $\mathbf{B}(\mathbf{x}) = \mathbf{B}(B_f,B_s;\mathbf{x})$ satisfies \eqref{eq-pcti}. (The elliptic case is obtained by setting $\phi = 0$.) The transmission form of \eqref{eq-parabolic} highlights the heterogeneity, %% \begin{subequations} \label{eq-par-transmission} \begin{eqnarray} \phi_f \dtt{u_f} - \nabla \cdot ( B_f \nabla u_{f} ) &=& f_f, \;\; \mathbf{x} \in \Omega_{f} \label{eq-par-multi-f} \\ \phi_s \dtt{u_i} - \nabla \cdot ( B_s \nabla u_{i} ) &=& f_i, \;\; \mathbf{x} \in \Omega_{is}, \;\; i=1,\ldots N_{\rm incl} \label{eq-par-multi-s} \\ u_i &=& u_f, \;\; \mathbf{y} \in \Gamma_i \label{eq-par-multi-bc} \\ (B_f \nabla u_f) \cdot \mathbf{n} &=& (B_s \nabla u_i ) \cdot \mathbf{n}, \;\; \mathbf{y} \in \Gamma_i \label{eq-par-multi-flux} %% \\ \label{eq-par-multi-extbc} u_f &=& 0, \;\; \mathbf{x} \in \partial \Omega, \\ %% \label{eq-par-multi-initf} u_f(\mathbf{x},0) &=& u_{f0}(\mathbf{x}), \;\; \mathbf{x} \in \Omega_f \\ \label{eq-par-multi-inits} u_i(\mathbf{x},0) &=& u_{i0}(\mathbf{x}), \;\; \mathbf{x} \in \Omega_{is}\,. \end{eqnarray} %% \end{subequations} %% Two issues of scale appear in the behavior of $u(\mathbf{x},t)$. The first multiscale feature is associated with the local spatial variations of $u(\mathbf{x},t)$ due to variations of $\mathbf{B}$. The second multiscale character is associated with the time scale of getting to stationary equilibrium determined by the proportions of $\phi,\mathbf{B}$ in $\Omega_f$ and $\Omega_s$ \cite{BP98,Panfilov}. Three situations of upscaling are described in the following. %%%% \subsection{Classical upscaling} \label{sec-classical} When the geometry and coefficients are assumed periodic, and the ratio of coefficients is independent of $\epsilon$, problem \eqref{eq-parabolic} is upscaled by homogenization \cite{BLP,SP,Allaire92,Horn97} to obtain an effective equation \mpcite{SP}{II.5.(5.13)}. This equation is satisfied by the limit $\tilde{u}$ of the local averages of solutions to \eqref{eq-parabolic} as the number of inclusions $N_{\rm incl} \to \infty$, or, equivalently, as the size of the inclusions $\epsilon \to 0$, and it has the same form, % \begin{eqnarray} \label{eq-upscale-parabolic} \tilde{\phi} \frac{\partial \tilde{u}}{\partial t} - \nabla \cdot (\mathbf{\tilde{B}} \nabla \tilde{u}) = \tilde{f}, \; \mathbf{x} \in \Omega, %\\ %\tilde{u}(\mathbf{x},0)=\tilde{u}_0(\mathbf{x}) \end{eqnarray} % but with the local averages $\tilde{f}(\mathbf{x}) = \lave{f}{\Omega_0}$ and {\em constant coefficients}. The first is the local average of $\phi$ defined by % \begin{eqnarray} \label{eq-def-tphi} \tilde{\phi} \equiv \lave{\phi}{\Omega_0} = \frac{1}{\abs{\Omega_0}} ( \phi_f \abs{\Omega_{0f}}+ \phi_s \abs{\Omega_{0s}}) = \phi_f \theta_f + \phi_s \theta_s\,. \end{eqnarray} % The effective constant matrix $\mathbf{\tilde{B}}=\mathbf{\tilde{B}}(B_f,B_s)$ is defined as \mpcite{SP}{II.2.9} \begin{subequations} \label{eq-upscaleBB} % \begin{eqnarray} \label{eq-upscaleB} (\mathbf{\tilde{B}})_{jk} &=& \frac{1}{\snorm{\Omega _0}{}} \int_{\Omega _0} B_{jm} (\mathbf{y}) (\delta_{mk} + \partial_m \tilde{\omega}_k(\mathbf{y})) dA, \end{eqnarray} % where $\tilde{\omega}_j, j=1,2$ is a solution of the local periodic {\em cell problem} \mpcite{JKO}{Ch.1(1.35), 1.45a)} % \begin{eqnarray} \label{eq-omega} \left\{ \begin{array}{ll} -\nabla \cdot \mathbf{B} \nabla \tilde{\omega}_j(\mathbf{y}) &= \nabla \cdot (\mathbf{B} \mathbf{e}_j),\;\; \mathbf{y} \in \Omega _0 \\ \tilde{\omega}_j {\rm\ is\ } & \Omega _0-{\rm\ periodic} \end{array} \right. \end{eqnarray} \end{subequations} % It is understood that the condition {\em ``$\tilde{\omega}_j$ is $\Omega _0$-{periodic}''} in \eqref{eq-omega} constrains not just the values of the function $\tilde{\omega}_j$ on $\partial \Omega_0$, but also its normal flux $\mathbf{B} \nabla \tilde{\omega}_j \cdot \mathbf{n}$. The effective equation \eqref{eq-upscale-parabolic} describes very well this ``low contrast'' case, and the fine-scale variations of coefficients have been averaged out of the problem. %%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Obstacle problem} Consider the special case with $\phi_s=0,B_s=0,$ and $f_i =0$. This arises when the coefficient $\mathbf{B}$ in \eqref{eq-parabolic} is replaced by $\mathbf{{B}^0}(\mathbf{x}) \equiv \mathbf{B}(B_f,0;\mathbf{x})$. It is also known as {\em perforated domain} case, see \cite{Allaire92} or references therein, or {\em soft-inclusion} in material science, see \cite{CiorPaulin79}, \mpcite{JKO}{Section 3.1}. Denote by $\tilde{u}^0(\mathbf{x},t)$ the corresponding solution of the upscaled model, %% \begin{eqnarray} \label{eq-upscale-parabolic-obstacle} {\phi}^* \frac{\partial \tilde{u}^0 }{\partial t} - \nabla \cdot( \mathbf{{B}^*} \nabla \tilde{u}^0) = {f^*}, \end{eqnarray} % where ${f^*} = \lave{f}{\Omega_0} = \theta_f \lave{f}{\Omega_{0f}}$. Note that there is no storage in $\Omega_{0s}$ and the local cells are impermeable. % %Note $\tilde{u}^0 = \lave{u^0}{\Omega_{0f}}$. % The upscaled constant coefficient $\mathbf{{B}^*}$ is given as before by % \begin{eqnarray} \label{eq-upscaleBtz} (\mathbf{{B}^*})_{jk} &=& \frac{1}{\snorm{\Omega_0}{}} \int_{\Omega_{0f}} B_{jm} (\mathbf{y}) (\delta_{mk} + \partial_m \tilde{\omega}^0_k(\mathbf{y})) dA. \end{eqnarray} % where $\tilde{\omega}^0(\mathbf{x})$ is defined as the solution over $\Omega_f$ of the cell problem % % \begin{eqnarray} \label{eq-obstacle} \left\{ \begin{array}{ll} -\nabla \cdot \nabla \tilde{\omega}^0_j(\mathbf{y}) = 0, \;&\mathbf{y} \in \Omega_{0f} \\ \nabla \tilde{\omega}^0_j(\mathbf{y}) \cdot \mathbf{n} = -\mathbf{e}_j \cdot \mathbf{n}, \; \mathbf{y} \in \Gamma_{fs}\\ \tilde{\omega}^0_j {\rm\ is\ }&\Omega _0-{\rm\ periodic}. \end{array} \right. \end{eqnarray} % We note that in some references \cite{JKO} the function $\tilde{\omega}^0$ is extended to all of $\Omega_0$ and that \eqref{eq-obstacle} agrees with \eqref{eq-omega}. Let us briefly compare the solutions $\tilde{u}$ and $\tilde{u}^0$ when $f_i = 0$. The former describes the evolution in time of local averages over $\Omega_0$, the latter is concerned with the averages over $\Omega_{0f}$. The former, at least formally, captures the transient storage in both regions $\Omega_{0f}$ and $\Omega_{0s}$, while the latter is only considered with evolution of storage in $\Omega_{0f}$. Finally, $\tilde{u}^0$ is the formal limit of $\tilde{u}$ as $B_s \to 0$. %%%%%%%%%%%% \subsubsection{Double porosity models} \label{sec-double-porosity} Here we follow the presentations \mpcite{ADH}{(3.4)}, \cite{Horn97}. and \mpcite{Allaire92}{Section 4}. In the last reference, this is called the {\em highly heterogeneous} case. Assume $f_i = 0$ and that the coefficient $B_s$ is scaled by $\epsilon^2$, {\em i.e.}, $B_s$ is replaced by $\epsilon^2 B_s$ throughout \eqref{eq-par-transmission}. This has the effect of balancing the decreasing size of the inclusions with the permeability to retain the coupling of the two regimes. (By contrast, in the obstacle case, the inclusions are impermeable.) The upscaled solution $u^*$ satisfies the equation % \begin{subequations} \label{eq-upscale-parabolic-ADH} \begin{eqnarray} \label{eq-upscale-parabolic-global} {\phi}^* \dtt{u^*} + \frac{1}{\abs{\Omega_0}} \int_{\Omega_{0s}} \phi_s \frac{\partial {u_s^*}_0 }{\partial t} dA - \nabla \cdot (\mathbf{{B}^*} \nabla u^*) = {f^*}, \; \mathbf{x} \in \Omega\,, \end{eqnarray} % in which the second term contains the solution ${u_s^*}_0$ of a problem on $\Omega_{0s}(\mathbf{x})$ at every $\mathbf{x} \in \Omega$, % \begin{eqnarray} \label{eq-upscale-parabolic-local} \phi_s \frac{\partial {u_s^*}_0}{\partial t} - \nabla \cdot (\mathbf{B_{s}} \nabla {u_s^*}_0) = 0,\; \mathbf{y} \in \Omega_{0s}(\mathbf{x}), \\ {u_s^*}_0 \vert_{\Gamma_0} = u^*(\mathbf{x},t). \end{eqnarray} % \end{subequations} % Here we make a distinction between global variable $\mathbf{x} \in \Omega$ and a local variable $\mathbf{y} \in \Omega_{0s}(\mathbf{x})$, where the local inclusion $\Omega_{0s}(\mathbf{x})$ is centered at $\mathbf{x}$. This is the double-porosity model for single-phase slightly compressible flow first developed in \cite{DougArb} and derived formally in \cite{ADH} and discussed in this form under the name ``distributed micro-structure model'' in \cite{ShowWalk91a,ShowWalk93}. The asymptotic expansions and the special weak convergence proof which preceded the two-scale convergence ideas appeared in \cite{DougArb,ADH,Arb89}. We emphasize that the constant tensor $\mathbf{{B}^*}$ which made the inclusions impermeable in the limit for the obstacle problem is precisely that obtained by the use of $\epsilon^2-$ scaling \cite{DougArb,ADH}. This scaling is important for the parabolic problem. Note that for the corresponding elliptic system discussed in \mpcite{Horn97}{pp18-22,Section 6.3.2,145}, a closer look at the problem reveals that ${u_s^*}_0\vert_{\Omega_{0s}} \equiv const$ and therefore the second term of \eqref{eq-upscale-parabolic-global} vanishes and the system is uncoupled. See Corollary~\ref{cor-ell} for further details and consequences. The $\epsilon^2$-scaling was first proposed in \cite{Vogt}; it can be also justified heuristically as in \cite{DougArb}. % On the contrary, in the model of molecular diffusion as in \cite{spagnuolo01}, no scaling is applied. A general family of scalings using $\epsilon^{\gamma}$ with $\gamma \neq 2$, was considered in \cite{AGP05}. See also \cite{PeterBohm05}, \cite{PanPiatRyb03}. In multiphase flow models in \cite{Bourgeat,Bourgeat97} the scaling by $\gamma <2$ was used while $\gamma=2$ was used in \cite{Arbogast}. Other models known as dual (double) porosity models similar to \eqref{eq-upscale-parabolic-ADH} had been proposed in applications \cite{BZK60,WarRoo63} prior to the introduction of formal homogenization methods and prior to \cite{DougArb}. These models have a structure similar to \eqref{eq-upscale-parabolic-ADH} but feature a general exchange/memory term $q^*(t)$ corresponding to a variety of problems on $\Omega_s$, which in turn are coupled by interface conditions. They all have the form % \begin{subequations} \label{eq-upscaledS-parabolic} \begin{eqnarray} \label{eq-upscaledS-global} {\phi}^* \frac{\partial u^*}{\partial t} + q^*- \nabla \cdot (\mathbf{{B}^*} \nabla u^*) = 0, \; \mathbf{x} \in \Omega\,. \end{eqnarray} % In some cases, these memory terms can be written as nonlocal Volterra terms in convolution form, % \begin{eqnarray} \label{eq-upscaledS-qterm} {q^*} = \tau \ast L (u^*)\,, \end{eqnarray} % \end{subequations} % where $L(u^*)$ is some differential operator and $\tau$ is some convolution kernel. In the double-porosity model \eqref{eq-upscale-parabolic-ADH}, the convolution kernel $\tau$ behaves, for small times, like $\tau(t) \approx \frac{1}{\sqrt{t}}$ \cite{HornSho90,P92} and acts on $L(u^*) \approx \dtt{u^*}$. In applications to flow and transport, it is the term ${q^*}$ that gives rise to the {\em tailing} in diffusion/dispersion \cite{Hagg95}. The dynamic and possibly delayed response of the cell $\Omega_{0s}$ and its effect on the global solution are quantified by representing it as a nonlocal-in time term \eqref{eq-upscaledS-qterm}. Depending on the singularity of $\tau$ at $0$ and on its long-term behavior one has more or less significant memory effects: for more singular kernels the memory is short-term, for less singular kernels the memory is long-term. In some cases considered in applications \eqref{eq-upscaledS-qterm} is shown to be multi-rate or its Prony series can be truncated \cite{P96} to capture effects most significant for the given case study \cite{Hagg95,Hagg04}. In general, a stochastic representation of \eqref{eq-upscaledS-qterm} has to be used to account for uncertainty. In \cite{spagnuolo01} the term $q^*(t) \approx \phi_s \frac{\partial u^*}{\partial t}$ when $\mathbf{D_s}$ represents molecular diffusion only; this amounts to identifying the convolution kernel $\tau$ with Dirac distribution $\delta_0$. Finally, it was shown recently that the different terms ${q^*}$ relate to a particular choice of $\gamma$ in $\epsilon$-scaling, see for example how the model in \cite{BZK60} was derived in \cite{ShowVis04}. In summary, the choice of which scaling and family of effective models to use should depend on both the ratio $\frac{B_f}{B_s}$ and the size $\epsilon_0$. In particular, since the limiting model \eqref{eq-upscale-parabolic-ADH} is derived assuming $\epsilon \to 0$, its use may be of limited value when $\epsilon_0$ is not small or $\frac{B_f}{B_s}$ is not large. Finally, the method of upscaling chosen for the elliptic flow equation \eqref{eq-flow} directly influences and sometimes inappropriately eliminates the advection terms in the upscaled transport equation \eqref{eq-transport}. These modeling issues motivated by experimental results in \cite{Hagg04} are combined in the following in order to describe a broad range of not well separated scales. \section{Discrete double porosity models with local approximations and projections} \label{sec-discrete} Here we construct a discrete version of \eqref{eq-upscaledS-parabolic}. Our development does not require taking the limit as $\epsilon\to 0$. The upscaled system contains as a special case discrete counterparts of the double-porosity parabolic model \eqref{eq-upscale-parabolic-ADH}. We begin with the exact discrete model \eqref{eq-par-transmission} and construct an upscaled model which uses a local approximation on the interfaces $\Gamma_s$; this refines the approximations in \cite{Arb89,Arb89r}. In order to conserve mass and, equivalently, to prevent creation of sources and sinks in the upscaled model, we derive a compatibility condition between the two approximations used in interface conditions \eqref{eq-par-multi-bc} and \eqref{eq-par-multi-flux}. The variational framework is used to obtain a consistent discrete form of ${q^*}$ in \eqref{eq-upscaledS-qterm} and to establish well-posedness. At the end in Section~\ref{sec-discrete-elliptic} we consider a stationary elliptic limit to the parabolic problem \eqref{eq-upscaledS-parabolic} and its discrete double-porosity counterpart. The result will be used in Section~\ref{sec-newmodel} for the flow part of the coupled flow-transport model \eqref{eq-flow}--\eqref{eq-transport}. \subsection{Variational form of exact discrete model} \label{sec-var-exact} The well-posedness of \eqref{eq-parabolic} or equivalently \eqref{eq-par-transmission} is standard; for double porosity models it was considered in a slightly different setup in \cite{Arb89,ShowWalk91a} and many other works. We begin by recalling the variational formulation which shows the dynamics are governed by an analytic semigroup. We use standard notation for Lebesgue spaces, $L^2(D)$, and Sobolev spaces, $H^s(D), H^r(\partial D),\; s,r \in \mathbb{R}$, and a standard symbol $\langle x',x\rangle \equiv \langle x',x\rangle_X$ to denote duality pairing between elements of a space $X$ and its dual $X'$ \cite{Adams}. We begin with the source-free case, $f = 0$. The weak formulation of the initial-boundary-value problem \eqref{eq-parabolic} is %% \begin{subequations} \label{eq-var-model} \begin{eqnarray} u(t) \in V: \quad \ddt{} \left(\phi\, u(t),w \right)_{H} + (\mathbf{B} \nabla u(t), \nabla w)_H &=&0, \;\;\forall w \in V, \\ u(\cdot,0) &=& u_0(\cdot) \in H \end{eqnarray} \end{subequations} %% where we have denoted by $(\cdot,\cdot)_H$ the scalar product in $H = L^2(\Omega)$ and set $V=H_0^1(\Omega)$. The operator on $H$ associated with the bilinear elliptic form $(\mathbf{B} \nabla u, \nabla w)_H$ on $V$ is known to be (the negative of) the generator of an analytic semigroup on $H$, and this leads to the following standard result. \begin{proposition} \label{prop-well} Let $u_0 \in H$. Then the initial-boundary-value problem \eqref{eq-var-model} has a unique solution $u(\cdot) \in C([0,\infty),H) \cap C^{\infty}((0,\infty),V)$. \end{proposition} The corresponding transmission form \eqref{eq-par-transmission} suggests an equivalent formulation which connects more naturally with the upscaled equation \eqref{eq-upscaledS-parabolic}. Note that the space $L^2(\Omega)$ can be identified with the product $L^2(\Omega_f) \times \prod_{i=1}^{N_{\rm incl}} L^2(\Omega_{is})$, and the Sobolev space $H_0^1(\Omega)$ is similarly identified with % \begin{multline} \{(w_f,\ivec{w_i}) \in H^1(\Omega_f) \times \prod_{i=1}^{N_{\rm incl}} H^1(\Omega_{is}): \phantom{satisfied by\ } \\ \eqref{eq-par-multi-bc} {\rm \ are\ satisfied\ by\ } w_f,w_i;\; \eqref{eq-par-multi-extbc} {\rm \ is\ satisfied\ by\ } w_f\}. \label{eq-def-v} \end{multline} % To get the weak form directly, we multiply \eqref{eq-par-multi-f} and \eqref{eq-par-multi-s} by the appropriate components of test functions, integrate over $\Omega_f$ and $\bigcup_i \Omega_{is}$ and add the equations. The weak form follows: find $\mathbf{U} =(u_f,\ivec{u_i}) \in V$ such that % \begin{multline} \label{eq-par-var} (\phi_f \dtt{u_f},v_f)_{L^2(\Omega_f)}+ \sum_i ((\phi_s \dtt{u_i},v_i)_{L^2(\Omega_{is})} \\ + (B_f \nabla u_{f}, \nabla v_f)_{L^2(\Omega_f)}+ \sum_i (B_s \nabla u_{i}, \nabla v_i)_{L^2(\Omega_{is})} \\ = \sum_i \int_{\Gamma_i} (B_f \nabla u_{f}v_f -B_s \nabla u_{i}v_i) \cdot \mathbf{n} ds,\;\; \forall \mathbf{V}=(v_f,\ivec{v_i}) \in V \end{multline} % Using \eqref{eq-def-v} we see that \eqref{eq-par-multi-flux} holds exactly when % \begin{eqnarray} \label{eq-ex-compatibility} \sum_i \int_{\Gamma_i} (B_f \nabla u_{f}v_f -B_s \nabla u_{i}v_i) \cdot \mathbf{n}=0,\;\; \forall \mathbf{V}=(v_f,\ivec{v_i}) \in V, \end{eqnarray} % that is, the right side of \eqref{eq-par-var} vanishes. The system takes the form % \begin{eqnarray} \label{eq-wk} \mathbf{U}(t) \in V: \quad (\phi \dtt{\mathbf{U}},\mathbf{V})_{H} + \mathcal{B}(\mathbf{U},\mathbf{V}) =0,\; \forall \mathbf{V} \in V, \end{eqnarray} % where $\mathcal{B}(\cdot,\cdot)$ is the continuous and coercive bilinear form % \[ \mathcal{B}(\mathbf{U},\mathbf{V}) \equiv (B_f \nabla u_{f}, \nabla v_f)_{L^2(\Omega_f)}+ \sum_i (B_s \nabla u_{i}, \nabla v_i)_{L^2(\Omega_{is})} \,, \quad \mathbf{U},\ \mathbf{V} \in V. \] % This is the essential ingredient for the proof of Proposition~\ref{prop-well}. \begin{remark} \rm One can add advection terms in the preceding. For a $\mathbf{b} \in (L^{\infty}(\Omega))^d$, the conclusion of Proposition~\ref{prop-well} holds \cite{Show-monotone} if the bilinear form is supplemented with first-order terms %% \begin{eqnarray} (\mathbf{B} \nabla u(t) + \mathbf{b} u(t), \nabla w)_H\,. \end{eqnarray} %% \end{remark} A condition analogous to \eqref{eq-ex-compatibility} is crucial in the derivation and analysis of the upscaled problems: it ensures that no internal sinks or sources are created in the process of upscaling. %%%%%%%%%% \subsection{Discrete upscaled model} \label{sec-discrete-upscaled} Now we propose and analyze an upscaled version of \eqref{eq-par-transmission} leading to a discrete analogue of \eqref{eq-upscaledS-parabolic}. The model includes as special cases the situations in \cite{Arb89,ShowWalk91a}, some elements of \cite{Arb89r,ClarkShow94}, the applications models \cite{Hagg95,Hagg00,wood.cmg}, and a discrete version of \eqref{eq-upscale-parabolic-ADH}. The gist of the construction is to identify properly the two-way coupling between the global equation \eqref{eq-upscaledS-global} and the local equation \eqref{eq-upscaledS-qterm}, and to define ${q^*}$ accordingly. Comparing \eqref{eq-par-transmission}, \eqref{eq-upscale-parabolic-obstacle}, and \eqref{eq-upscale-parabolic-ADH}, we expect an upscaled form %% \begin{subequations} \label{eq-par-model} \begin{eqnarray} {\phi}^* \frac{\partial u^*}{\partial t} + {q^*}(\mathbf{x},t) - \nabla \cdot ( \mathbf{{B}^*} \nabla u^* ) &=& 0, \;\; \mathbf{x} \in \Omega, \label{eq-par-model-global-pde} \\ \phi_{s} \frac{\partial {u_s^*}_{i}}{\partial t} - \nabla \cdot ( B_s \nabla {u_s^*}_{i} ) &=& 0, \;\; \mathbf{y} \in \Omega_{is} \label{eq-par-model-block} \\ {u_s^*}_i \vert_{\Gamma_i} &=& (\mathbf{\Pi} u^*)_i \label{eq-par-model-bc} \\ \label{eq-par-model-extbc} u^* &=& 0, \;\; \mathbf{x} \in \partial \Omega, \\ \label{eq-par-model-initf} u^* (\mathbf{x},0) &=& u^{0}(\mathbf{x}), \;\; \mathbf{x} \in \Omega \\ \label{eq-par-model-inits} {u_s^*}_i(\mathbf{y},0) &=& u_i^{0}(\mathbf{y}), \;\; \mathbf{y} \in \Omega_{is} \end{eqnarray} % \end{subequations} %% where $\mathbf{\Pi}$ and ${q^*}(\mathbf{x},t)$ are to be defined below. We intend for \eqref{eq-par-model-bc} to approximate \eqref{eq-par-multi-bc} while \eqref{eq-par-multi-flux} is realized by the flux term ${q^*}$ in \eqref{eq-par-model-global-pde}. Clearly one could come up with many heuristic reasonable {\em independent} choices for $\mathbf{\Pi}$ and ${q^*}$. We argue that in order to conserve mass and preserve the variational structure exhibited for the exact problem~\eqref{eq-par-transmission}, these choices cannot be made independently. We show below that ${q^*}$ must contain the dual operator to $\mathbf{\Pi}$. We seek the weak formulation of \eqref{eq-par-model} which will lead to the appropriate compatibility condition. The solution to \eqref{eq-par-model} is an $(N_{\rm incl} +1)$-vector % \[{\mathbf{U}^*}\equiv (u^*,({u_s^*}_i)_i) \in {\mathcal \mathbf{H}} \equiv L^2(\Omega) \times \prod_{i=1}^{N_{\rm incl}} L^2(\Omega_{is}).\] % Note that ${\mathcal \mathbf{H}}$ is much more than and cannot be identified with $H = L^2(\Omega)$. The weak solution ${\mathbf{U}^*}$ should also belong to $X \times \mathbf{Y}$ with % \[X=H_0^1(\Omega), \; \mathbf{Y}= \prod_{i=1}^{N_{\rm incl}} H^1(\Omega_{is}).\] % The operator $\mathbf{\Pi}$ which appears in \eqref{eq-par-model-bc} takes values in $\Gamma_i$, so for $w \in X$ its $i$-th component $(\mathbf{\Pi} w)_i \in Z_i \equiv H^{1/2}(\Gamma_i)$, the restrictions (traces) to $\Gamma_i$ of functions from $\mathbf{Y}$. Thus we define $ \mathbf{\Pi}: X \rightarrow \mathbf{Z}, $ where \[\mathbf{Z} \equiv \prod_{i=1}^{N_{\rm incl}} Z_i = \prod_{i=1}^{N_{\rm incl}} H^{1/2}(\Gamma_i).\] % Recall also the characterization of the dual spaces, \[X' = H^{-1}(\Omega), \;\;\; \mathbf{Z}'=\prod_{i=1}^{N_{\rm incl}} H^{-1/2}(\Gamma_i).\] We embed the interface conditions \eqref{eq-par-model-bc} in the definition of the energy space ${\mathbf{V}}$, \[ {\mathbf{V}}\equiv \{(w,\ivec{w_i}) \in X \times \mathbf{Y}: % w_i \vert_{\Gamma_i} = (\mathbf{\Pi} w)_i, \forall i\}, \] where the constraint is understood in the sense of traces. % For a sufficiently regular solution ${\mathbf{U}^*}=(u^*,\ivec{{u_s^*}_i}) \in {\mathbf{V}}$, the weak form of \eqref{eq-par-model} is obtained by multiplying \eqref{eq-par-model-global-pde}, \eqref{eq-par-model-block} by the appropriate components of a test function $\mathbf{W}=(w,\ivec{w_i})\in {\mathbf{V}} $ and integrating over $\Omega$ and $\Omega_{is}$, respectively, % \begin{eqnarray} % \int_{\Omega} {\phi}^* \frac{\partial u^*}{\partial t} w (\mathbf{x}) dA + \int_{\Omega} {q^*}(\mathbf{x},t) w(\mathbf{x}) dA +\int_{\Omega} \mathbf{{B}^*} \nabla u^* \nabla w(\mathbf{x}) dA =0\,, \label{eq-g} \\ \int_{\Omega_{is}} \phi_{s} \frac{\partial {u_s^*}_i}{\partial t} (\mathbf{y},t) w_i(\mathbf{y}) dA + \int_{\Omega_{is}} B_s \nabla {u_s^*}_i (\mathbf{y},t)\nabla w_i(\mathbf{y})dA \nonumber \\ = \int_{\Gamma_i} q_i(\mathbf{s},t) w_i(\mathbf{s}) dS\,, \label{eq-i} \end{eqnarray} % where we have denoted the flux $q_i \in Z_i'$ by % $ q_i (\mathbf{s},t) \equiv (B_s \nabla {u_s^*}_i \cdot \mathbf{n})(\mathbf{s}),\; \mathbf{s} \in \Gamma_i $. % Sum \eqref{eq-i} over $i$, add the result to \eqref{eq-g}, and use the relation $w_i=(\mathbf{\Pi} w)_i$ which was embedded in the definition of space ${\mathbf{V}}$. The resulting system is %% \begin{multline} \ddt{} \left(\phi {\mathbf{U}^*},\mathbf{W} \right)_{\mathbf{H}} + \mathcal{B}({\mathbf{U}^*},\mathbf{W}) = \\ \left[ \sum_i \int_{\Gamma_i} q_i(\mathbf{s},t) (\Pi w )_i(\mathbf{s}) dS - \int_{\Omega} {q^*}(\mathbf{x},t) w(\mathbf{x}) dA\right], \;\;\forall \mathbf{W} \in \mathbf{V}, \label{eq-model-wp} \end{multline} %% where $\mathcal{B}$ is the positive definite bilinear form % \begin{eqnarray} \mathcal{B}({\mathbf{U}^*},\mathbf{W}) \equiv (\mathbf{{B}^*} \nabla u^*, \nabla w)_{\Omega} + \sum_i (B_s \nabla {u_s^*}_{i}, \nabla w_i)_{\Omega_{is}}. \label{p-d-f} \end{eqnarray} % If we have %% \begin{eqnarray} \int_{\Omega} {q^*}(\mathbf{x},t) w(\mathbf{x}) dA = \sum_i \int_{\Gamma_i} q_i(\mathbf{s},t) (\Pi w )(\mathbf{s}) dS, \quad \mathbf{W} \in \mathbf{V}, \label{eq-dual0} \end{eqnarray} %% then the weak formulation of \eqref{eq-par-model} is \eqref{eq-wk} with \eqref{p-d-f}. The condition \eqref{eq-dual0} can be interpreted as a statement ensuring conservation of mass. It requires that the term ${q^*}$ be compatible with the collection of fluxes $\mathbf{q}=\ivec{q_i}$ out of $\Omega_s$ so that the right side of \eqref{eq-model-wp} vanishes and there are no spurious sources or sinks in the model. In the general case, since $\ivec{u_i} \in \mathbf{Y}$, the fluxes $\ivec{q_i}$ across $\partial \Omega_{is} \equiv \Gamma_i$ are in $\mathbf{Z}'$. On the other hand, $\ivec{w_i\vert_{\Gamma_i}} \in \mathbf{Z}$, so $\sum_{i=1}^{N_{\rm incl}} \int_{\Gamma_i} q_i w_i$ is understood as the duality pairing between $\mathbf{Z}$ and $\mathbf{Z}'$, $\langle \mathbf{q}, \mathbf{w}\rangle_{\mathbf{Z}}$. Similarly, the left side of \eqref{eq-dual0} for $w \in X$ and ${q^*} \in X'$ is % $\langle {q^*},w \rangle_X$, so we obtain condition \eqref{eq-dual0} in the form % \begin{eqnarray} \label{eq-dual1} \langle {q^*} ,w \rangle_X = \langle \mathbf{q}, \ivec{(\mathbf{\Pi} w)_i} \rangle_{\mathbf{Z}} = \langle \mathbf{q}, \mathbf{\Pi} w \rangle_{\mathbf{Z}}. \end{eqnarray} %% This shows that % \begin{eqnarray} {q^*} = \Pi' \mathbf{q}, \label{eq-dual} \end{eqnarray} % where $\Pi': \mathbf{Z}' \rightarrow X'$ is the operator dual to $\mathbf{\Pi}: X \rightarrow \mathbf{Z}$, % \begin{eqnarray} \langle \Pi' \mathbf{q} ,w \rangle_X \equiv \langle \mathbf{q}, \mathbf{\Pi} w \rangle_{\mathbf{Z}}, \;\; \mathbf{q} \in \mathbf{Z}', w \in X. \label{eq-dual-operator} \end{eqnarray} % With this characterization, the model \eqref{eq-par-model} is rewritten precisely as %% \begin{subequations} \label{eq-par-modelf} % extended \begin{eqnarray} {\phi}^* \frac{\partial u^*}{\partial t} + \Pi' \left( \ivec{B_s \nabla {u_s^*}_i \cdot \mathbf{n}_i} \right) - \nabla \cdot ( \mathbf{{B}^*} \nabla u^* ) &=& 0, \;\; \mathbf{x} \in \Omega, \label{eq-par-modelf-global-pde} \\ \phi_{s} \frac{\partial {u_s^*}_i}{\partial t} - \nabla \cdot ( B_s \nabla {u_s^*}_i ) &=& 0, \;\; \mathbf{y} \in \Omega_{is} \label{eq-par-modelf-block} \\ {u_s^*}_i \vert_{\Gamma_i} &=& (\mathbf{\Pi} u^*)_i\,, \label{eq-par-modelf-bc} \end{eqnarray} % \begin{eqnarray} \label{eq-par-modelf-extbc} u^* (\mathbf{x},t) &=& 0, \;\; \mathbf{x} \in \partial \Omega, \\ \label{eq-par-modelf-initf} u^* (\mathbf{x},0) &=& u^{0}(\mathbf{x}), \;\; \mathbf{x} \in \Omega, \\ \label{eq-par-modelf-inits} {u_s^*}_i(\mathbf{y},0) &=& u_i^{0}(\mathbf{y}), \;\; \mathbf{y} \in \Omega_{is}. \end{eqnarray} %% \end{subequations} % We summarize the preceding in the following. %% \begin{proposition} \label{prop-var-pi} The weak formulation of the upscaled discrete system \eqref{eq-par-model} is \eqref{eq-par-modelf}, and it is a well-posed initial-boundary-value problem. The evolution of the solutions to the model \eqref{eq-par-modelf} is governed by an analytic semigroup on $\mathbf{H}$. \end{proposition} \subsection{The Operators $\mathbf{\Pi}$} \label{sec-pi} It remains to define $\mathbf{\Pi}$ in the upscaled model \eqref{eq-par-modelf}. Its role is to provide an appropriate approximation to the boundary values in \eqref{eq-par-modelf-bc} for which the resulting discrete upscaled model is accurate and numerically tractable. We consider only at most linear polynomial approximations; for $m=0,1$, we use the space $P_m$ of polynomials of degree at most $m$. Denote by $Z_i^m = P_m(\Gamma_i)$ these polynomials regarded as functions on $\Gamma_i$, and corresponding subspaces of $\mathbf{Z}$ for the operator domain, %% \begin{eqnarray} {\mathbf{Z}}^m \equiv \prod_{i=1}^{N_{\rm incl}} Z_i^m, \quad \label{eq-map} \mathbf{\Pi}_m:X \rightarrow \mathbf{Z}^m. \end{eqnarray} % Piecewise constant and affine approximations were used in \cite{Arb89,ShowWalk91a} and in \cite{Arb89r,ClarkShow94,spagnuolo01}, respectively. The discussion of corresponding maps $\mathbf{\Pi}_m:X\rightarrow \mathbf{Z}^m$ for $m=0$ and $m=1$ are carried out in Section~\ref{sec-constant} and \ref{sec-affine}, respectively. The calculations of $\mathbf{\Pi}_m'$ lead to moments of $q_i \in Z_i'$. The zero'th and first order moments are %% \begin{eqnarray} \label{eq-moment0} M_i^0 (q_i) \equiv \frac{1}{\abs{\Omega_i}} \langle q_i,1\rangle _{Z_i}, \;\; \mathbf{M}_i^1(q_i)&\equiv & \frac{1}{\snorm{\Omega_i}{}} \langle q_i, (\mathbf{s} - \mathbf{x}_i^C)\rangle _{Z_i}, \;\; q_i \in Z_i'. \end{eqnarray} %% For smoother $q_i$ they are given by %% \begin{eqnarray} \nonumber M_i^0 (q_i) = \frac{1}{\abs{\Omega_i}} \int_{\Gamma_i} q_i(\mathbf{s}) dS, \;\; \mathbf{M}_i^1(q_i)&=& \frac{1}{\snorm{\Omega_i}{}} \int_{\Gamma_i} q_i(s) (\mathbf{s} - \mathbf{x}_i^C) dS, \;\; q_i \in L^2(\Gamma_i). \end{eqnarray} %% Recall that we have identified a constant or linear polynomial as a function on $\Gamma_i$, hence, an element of $Z_i$, so these definitions make sense. We will also use the following consequence of the Green's formula applicable to the first moments. % \begin{lemma} \label{lem-Green} For any smooth region $D$, smooth $\mathbf{v}=(v_1,v_2)$ and the centroid $\mathbf{x}^C_D \in D$ of $D$, we have, for $k=1,2$ \[\int_{D} (\nabla \cdot \mathbf{v}) (x_k-(\mathbf{x}^C_D)_k) dA = \int_{\partial D} \mathbf{v} \cdot \mathbf{n} (x_k-(\mathbf{x}^C_D)_k)dS - \int_D v_k dA. \] \end{lemma} \subsubsection{Piecewise constant approximations $\Pi_0$} \label{sec-constant} Here we discuss the local constant approximations. Define $\mathbf{\Pi}_0: X \rightarrow \mathbf{Z}^0$ as the local averages over $\Omega_i$, % \begin{eqnarray} \mathbf{\Pi}_0 w &\equiv& \ivec{(\mathbf{\Pi}_0 w)_i}, \\ (\mathbf{\Pi}_0 w)_i(\mathbf{s}) &\equiv& {\lave{w}{i}} = \frac{1}{\abs{\Omega_i}} \int_{\Omega_i} w(\mathbf{x}) dA\,, \;\; \mathbf{s} \in \Gamma_i. \end{eqnarray} %% The dual operator $\Pi_0' : (\mathbf{Z}^0)' \rightarrow X'$ is computed as follows, assuming $\mathbf{q}$ is sufficiently smooth, as it will be in the system. We have %% \begin{multline*} \langle \mathbf{q}, \mathbf{\Pi}_0 w\rangle_{\mathbf{Z}} = \sum_i \int_{\Gamma_i} (\Pi_0 w)_i (s) q_i(s) dS = \sum_i \int_{\Gamma_i} \frac{1}{\abs{\Omega_i}} \left(\int_{\Omega} \hat{\chi}_i(\mathbf{x}) w(\mathbf{x}) dA \right) q_i(s) dS \\ = \int_{\Omega} w(x) \left( \sum_i \hat{\chi}_i(\mathbf{x}) \frac{1}{\abs{\Omega_i}} \int_{\Gamma_i} q_i(s) dS \right) dA = \int_{\Omega} w(x) \left( \sum_i \hat{\chi}_i(\mathbf{x}) M_i^0 (q_i) \right) dA \end{multline*} %% so we obtain the following characterization of $(\Pi_0' \mathbf{q})$; its second part follows by divergence theorem. % \begin{lemma} \label{lem-ppi0} For $\mathbf{q}$ with integrable components we identify pointwise, in the sense of distributions %% \begin{eqnarray} \label{eq-pi0} (\Pi_0' \mathbf{q})(\mathbf{x}) = \sum_i \hat{\chi}_i(\mathbf{x}) M_i^0 (q_i). \end{eqnarray} % Additionally, let $q_i = \mathbf{q}_i \cdot \mathbf{n}$ for some smooth $\mathbf{q}_i$. In this case %% \begin{eqnarray} \label{eq-moment0div} (\Pi_0' (\mathbf{q}_i\cdot \mathbf{n}))(\mathbf{x}) = \sum_i \hat{\chi}_i(x) \frac{1}{\abs{\Omega_i}} \int_{\Omega_{is}} \nabla \cdot \mathbf{q}_i(\mathbf{x}) dA. \end{eqnarray} \end{lemma} % Now we want to point out the formal connection between the discrete and continuous models obtained from the following simple considerations expressing the convergence of piecewise constant approximations. Let $\Omega$ be fixed; as $N_{\rm incl} \to \infty$, we have that $diam(\Omega_i) \to 0$. For $u \in L^2(\Omega)$, define $\mathcal{T}_{N_{\rm incl}} (u)(\mathbf{x}) \equiv \sum_{i=1}^{N_{\rm incl}} \hat{\chi}_i(\mathbf{x}) \frac{1}{\abs{\Omega_i}} \int_{\Omega_i} u(\mathbf{y}) dA$. Then %% \[ (\mathcal{T}_{N_{\rm incl}} (u)(\mathbf{x}) )^2 = \sum_{i=1}^{N_{\rm incl}} \hat{\chi}_i(\mathbf{x}) \frac{1}{\abs{\Omega_i}^2}\left( \int_{\Omega_i} u(\mathbf{y}) dA\right)^2 \leq \sum_{i=1}^{N_{\rm incl}} \hat{\chi}_i(\mathbf{x}) \frac{1}{\abs{\Omega_i}}\int_{\Omega_i} (u(\mathbf{y}))^2 dA, \] so $\norm{\mathcal{T}_{N_{\rm incl}} (u)}{L^2(\Omega)} \leq \norm{u}{L^2(\Omega)}$. For a smooth $u \in C^0(\Omega)$, we have \[ \abs{\mathcal{T}_{N_{\rm incl}} (u)(\mathbf{x}) - u(\mathbf{x}) } \leq %\max_{1 \leq i \leq N_{\rm incl}} %\frac{1}{\abs{\Omega_i}} \int_{\Omega_i} %\abs{u(\mathbf{x}) - u(\mathbf{y})} dA \leq \max_{1 \leq i \leq N_{\rm incl}, \mathbf{y},\mathbf{x} \in \Omega_i} \abs{u(\mathbf{x}) - u(\mathbf{y})} \] which converges to $0$. Thus $\lim_{N_{\rm incl} \to \infty} \mathcal{T}_{N_{\rm incl}} (u)=u \;\; {\rm\ in\ } L^2(\Omega)$, and by density of $C^0(\Omega)$ in $L^2(\Omega)$ we have %% \begin{eqnarray} \label{eq-limit} \lim_{N_{\rm incl} \to \infty} \sum_{i=1}^{N_{\rm incl}} \hat{\chi}_i(\cdot) (\mathbf{\Pi}_0 u)_i = u(\cdot) \;\; {\rm\ in\ } L^2(\Omega) \end{eqnarray} %% for each $u \in L^2(\Omega)$. Analogous properties hold for derivatives of a function in $H^1(\Omega)$. \begin{remark} \rm We recognize now that the double-porosity model \eqref{eq-upscale-parabolic-ADH} is the limit, in the sense of \eqref{eq-limit}, of the model \eqref{eq-par-modelf} with the choice of $\mathbf{\Pi}=\mathbf{\Pi}_0$. \end{remark} \subsubsection{Local affine projections $\mathbf{\Pi}_1$} \label{sec-affine} Now we consider local affine approximations associated with the operator $\mathbf{\Pi}_1 : H_0^1(\Omega) \rightarrow \mathbf{Z}^1$. These are needed to capture the effects of advection and secondary diffusion in upscaled coupled models; see Section~\ref{sec-newmodel}. There are many ways to define local affine approximations. One way is to use local Taylor approximations \cite{spagnuolo01}, but this requires extra smoothness. Another way proposed in \cite{Arb89r} which is mass conservative is to use local least-squares projections, see Remark~\ref{rem-lsq}. % Here we use affine approximations based on the moments defined by \eqref{eq-moment0}. We refer to these as local $H^1(\Omega_i)$-projections. % Define, for $i=1,2 \ldots N_{\rm incl}$, and $\mathbf{s} \in \Gamma_i$ %% \begin{eqnarray} \label{eq-pi-1} (\mathbf{\Pi}_1 w)_i(\mathbf{s}) &\equiv& (\mathbf{\Pi}_0 w)_i + (\mathbf{\Pi}_0 \nabla w)_i \cdot (\mathbf{s}-\mathbf{x}_i^C) \\ &= &\frac{1}{\snorm{\Omega_i}{}} \left( \int_{\Omega_i} w(\mathbf{y})\,dA + \sum_{k=1}^2 \left[\int_{\Omega_i} \partial_k w(\mathbf{y})\,dA\right]\, (s_k - (\mathbf{x}_i^C)_{k})\ \right),\; \mathbf{s} \in \Gamma_i. \nonumber \end{eqnarray} % %We notice that $\Pi_1$ is a projection since $\Pi_1 (\Pi_1 w) = \Pi_1 w$. To compute the dual $\Pi_{1}'$ to $\mathbf{\Pi}_1$ formally, let $\mathbf{q}$ have integrable components $q_i, \forall i =1,\ldots N_{\rm incl}$, and let $w \in X$. Then % \begin{multline} \label{eq-pidual} \langle {\Pi_1' \mathbf{q} },w \rangle_X = \langle \mathbf{q},{\mathbf{\Pi}_1 w} \rangle_{\mathbf{Z}} = \int_{\Omega} \sum_i \hat{\chi}_i(\mathbf{x}) \left(\frac{1}{\snorm{\Omega_i}{}} \int_{\Gamma_i} q_i(\mathbf{s}) dS\right) w(\mathbf{x}) dA \\ + \int_{\Omega} \sum_i \hat{\chi}_i(\mathbf{x}) \left(\frac{1}{\snorm{\Omega_i}{}} \int_{\Gamma_i} q_i(\mathbf{s}) (\mathbf{s} - \mathbf{x}^C) dS\right) \cdot \nabla w(\mathbf{x}) dA. \end{multline} %% We have thus shown the first part of the following Lemma. %% \begin{lemma} \label{lem-ppi1} Let $\mathbf{q}$ have integrable components $q_i,\forall i=1,\ldots N_{\rm incl}$. Then, in the sense of distributions in $H^{-1}(\Omega)$, the operator $\Pi_1'(\mathbf{q})(\mathbf{x})$ is characterized ``pointwise'' by % \begin{eqnarray} \label{eq-pi1} \Pi_1'(\mathbf{q})(\mathbf{x}) = \sum_i M_i^0 (q_i) {\hat{\chi}_i(\mathbf{x})} - \nabla \cdot \sum_i M_i^1 (q_i) {\hat{\chi}_i(\mathbf{x})}. \end{eqnarray} %% Note the second term in this equation is a collection of scaled line sources; indeed a member of $X'$. %% \\ Furthermore, let $q_i = \mathbf{q}_i \cdot \mathbf{n}$ for some $\mathbf{q}_i$ smooth on $\Omega_{is}$. Then %% \begin{multline} \label{eq-moment1div} \Pi_1'( (\mathbf{q}_i \cdot \mathbf{n})_i)(\mathbf{x}) = \sum_i \hat{\chi}_i(\mathbf{x}) \frac{1}{\abs{\Omega_i}} \int_{\Omega_{is}} \nabla \cdot \mathbf{q}_i(\mathbf{x}) dA \\ - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} \frac{1}{\abs{\Omega_i}} \int_{\Omega_{i,s}} (\nabla \cdot \mathbf{q}_i) (\mathbf{y} - \mathbf{x}^C) dA\, - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} \frac{1}{\abs{\Omega_i}} \int_{\Omega_{is}} \mathbf{q}_i dA\,. \end{multline} \end{lemma} %% \begin{proof} The proof of the second part follows easily if we apply Lemma~\ref{lem-Green} and calculate further from \eqref{eq-pi1} % \begin{eqnarray*} {\mathbf{M}_i^1(\mathbf{q}_i \cdot \mathbf{n})} = \frac{1}{\abs{\Omega_i}} \int_{\Omega_{is}} \left( (\nabla \cdot \mathbf{q}_i) (\mathbf{y} - \mathbf{x}^C) + \mathbf{q}_i \right) dA. \end{eqnarray*} % Next, incorporate \eqref{eq-moment0div} to get % \begin{multline*} \Pi_1'( (\mathbf{q}_i \cdot \mathbf{n})_i)(\mathbf{x}) = \sum_i M_i^0 (\mathbf{q}_i \cdot \mathbf{n}) {\hat{\chi}_i(\mathbf{x})} \\ - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} \frac{1}{\abs{\Omega_i}} \int_{\Omega_{is}} (\nabla \cdot \mathbf{q}_i) (\mathbf{y} - \mathbf{x}^C) dA\, - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} \frac{1}{\abs{\Omega_i}} \int_{\Omega_{is}} \mathbf{q}_i dA \end{multline*} %% from which \eqref{eq-moment1div} follows. % Note again that $\Pi_1'$ is a collection of line sources and thus requires test functions to be locally from $H^1$. \end{proof} \begin{corollary} \label{cor-par-pi1} The global equation of discrete upscaled parabolic model for \eqref{eq-par-model} constructed via compatibility condition \eqref{eq-dual} and given by \eqref{eq-par-modelf-global-pde}, with the choice $\mathbf{\Pi}=\mathbf{\Pi}_1$, takes the form %% \begin{multline} \label{eq-par-modelf-ppi1-global-pde} % {\phi}^* \frac{\partial u^*}{\partial t} + \sum_i \hat{\chi}_i(\mathbf{x}) \frac{1}{\abs{\Omega_i}} \int_{\Omega_{is}} \nabla \cdot (\mathbf{B_{s}} \nabla {u_s^*}_i) dA \\ - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} \frac{1}{\abs{\Omega_i}} \int_{\Omega_{i,s}} (\nabla \cdot (\mathbf{B_{s}} \nabla {u_s^*}_i) ) (\mathbf{y} - \mathbf{x}^C) dA\, \\ - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} \frac{1}{\abs{\Omega_i}} \int_{\Omega_{is}} (\mathbf{B_{s}} \nabla {u_s^*}_i) dA - \nabla \cdot ( \mathbf{{B}^*} \nabla u^* ) = 0, \;\; \mathbf{x} \in \Omega, \end{multline} %% where ${u_s^*}_i$ is the solution to \eqref{eq-par-modelf-block}--\eqref{eq-par-modelf-bc}. \end{corollary} In Section~\ref{sec-convolution} we provide calculations which partially decouple the global and local equations, and we show the continuous limit to \eqref{eq-par-modelf-ppi1-global-pde}. We close this section with a note on local least squares projections. % \begin{remark} \label{rem-lsq} \rm In some geometries the functions $(1,(\mathbf{x}-\mathbf{x}_i^C)_1,(\mathbf{x}-\mathbf{x}_i^C)_2)$ may be $L^2(\Omega_{is})$-orthogonal. Then, modulo proper normalization, $\mathbf{\Pi}_1$ is close to local $L^2(\Omega_{is})$-projections $\tilde{\mathbf{\Pi}}_1$ onto affines. The latter has a local orthonormal basis, preserves mass, and was used in a computational model in \cite{Arb89r}. \end{remark} For completeness we provide the explicit definition of $\tilde{\mathbf{\Pi}}_1$ and calculate $\tilde{\Pi}_1'$. Let $(\phi_i^0,\phi_i^1,\phi_i^2)$ be a local orthonormal basis functions spanning $Y_i^1$. We define $\tilde{\mathbf{\Pi}}_1$ componentwise as the $L^2(\Omega_i)$-projection onto affines $(\tilde{\mathbf{\Pi}}_1 w)_i(\mathbf{x}) \equiv \sum_k w_i^k \phi_i^k(\mathbf{x})$ where the coefficients $w_i^k$ are computed in a standard way via $w_i^k=\int_{\Omega_i} w(\mathbf{x}) \phi_i^k(\mathbf{x}) dA$. A calculation similar to the one for $\Pi_1'$ reveals that %% \begin{eqnarray*} (\tilde{\Pi}'_1 \mathbf{q})(\mathbf{x}) \equiv \sum_i \hat{\chi}_i(\mathbf{x}) \sum_k q_i^k \phi_i^k(\mathbf{x}). \end{eqnarray*} %%' where $q_i^k \equiv \int_{\Gamma_i} \phi_i^k(s) q_i(s) dS $. % It is apparent that $\tilde{\Pi}_1' \mathbf{q}$ is a globally discontinuous distribution of local polynomials. This is compatible with the fact that $\tilde{\mathbf{\Pi}}_1$ can be actually defined for (extended to) all $L^2(\Omega)$ functions not just those from $X$ and therefore its dual $\tilde{\Pi}'$ can be restricted to $L^2$ functions. The use of $\tilde{\mathbf{\Pi}}_1,\tilde{\Pi}_1'$ is straightforward. However, we were not able to find an interpretation of Fourier coefficients $q_i^k$ that would be useful in our subsequent model development. On the other hand, our $H^1(\Omega_i)$-projections while formally different from local least-squares projections may not be very different quantitatively. %%%%%%%%%%% \subsection{Discrete upscaled elliptic problem} \label{sec-discrete-elliptic} We close this Section with remarks on the elliptic counterpart of \eqref{eq-par-modelf} without local sources. % Consider the elliptic counterpart of \eqref{eq-par-transmission} with $f_i = 0$ and its discrete upscaled version. The well posedness of the variational formulation follows from the arguments above, Section~\ref{sec-discrete-upscaled}; see standard variational setting, e.g., \cite{ShowBook}. % Similar arguments as those for parabolic problems lead to a compatibility condition \eqref{eq-dual}. The discrete upscaled version for the elliptic case is obtained analogously as %% \begin{subequations} \label{eq-ell-modelf} % extended \begin{eqnarray} \Pi' \left( \ivec{B_s \nabla {u_s^*}_i \cdot \mathbf{n}_i} \right) &-& \nabla \cdot ( \mathbf{{B}^*} \nabla u^* ) = {f^*}, \;\; \mathbf{x} \in \Omega, \label{eq-ell-modelf-global-pde} \\ % - \nabla \cdot ( B_s \nabla {u_s^*}_i ) &=& 0, \;\; \mathbf{y} \in \Omega_{is},\; i=1,\ldots N_{\rm incl} \label{eq-ell-modelf-block} \\ % {u_s^*}_i \vert_{\Gamma_i} &=& (\mathbf{\Pi} u^*)_i. \label{eq-ell-modelf-bc} \\ \label{eq-ell-modelf-extbc} u^* (\mathbf{x}) &=& 0, \;\; \mathbf{x} \in \partial \Omega\,. \end{eqnarray} %% \end{subequations} We examine now the term ${q^*} \equiv \Pi' \left( \ivec{B_s \nabla {u_s^*}_i \cdot \mathbf{n}_i} \right)$ using Lemma~\ref{lem-ppi0} and Lemma~\ref{lem-ppi1}. The values in the boundary condition \eqref {eq-ell-modelf-bc} are either a constant $(\mathbf{\Pi}_0 u^*)_i$ or an affine function $(\mathbf{\Pi}_1 u^*)_i$. In either case the solution to \eqref{eq-ell-modelf-block} with \eqref{eq-ell-modelf-bc} is actually equal to the boundary data, constant or affine, respectively. Therefore, the flux of that solution $\mathbf{q}_i \equiv B_s \nabla {u_s^*}_i$ is either zero or constant, respectively. In the latter case, that constant equals $B_s (\mathbf{\Pi}_0 \nabla u^*)_i$. As a consequence, the zero'th moment of $\ivec{q_i}$ with $q_i \equiv \mathbf{q}_i \cdot \mathbf{n}$ which is (part of) the expression for $\Pi' \left( \ivec{q_i} \right)$ vanishes (see \eqref{eq-moment0div} and \eqref{eq-moment1div}) for both piecewise constant and piecewise affine cases. For the affine case we additionally find that the second term in \eqref{eq-moment1div} vanishes as well. For the third term in \eqref{eq-moment1div}, since $\mathbf{q}_i$ is locally constant, we get %% \begin{multline} \label{eq-third} - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} \frac{1}{\abs{\Omega_i}} \int_{\Omega_{is}} \mathbf{q}_i dA = - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} \frac{\abs{\Omega_{is}} }{\abs{\Omega_i}} B_s \nabla u^*_i \\ = - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} \theta_s B_s \nabla (\mathbf{\Pi}_0 \nabla u^* )_i\,. \end{multline} %% The consequences of these observations are summarized below. % \begin{corollary} \label{cor-ell} If a constant approximation associated with $\mathbf{\Pi}_0$ is used in the boundary condition \eqref{eq-ell-modelf-bc} for solutions to \eqref{eq-ell-modelf-block}, then (i) the flux $\mathbf{q}_i = B_s \nabla u^*_i$ equals $\mathbf{q}_i \equiv \mathbf{0}$, (ii) the source term ${q^*} \equiv 0$, and (iii) the model \eqref{eq-ell-modelf-global-pde} is just the obstacle problem \eqref{eq-upscale-parabolic-obstacle} with ${\phi}^* = 0$. % However, if the affine approximation $\mathbf{\Pi}_1$ is used in \eqref{eq-ell-modelf-bc}, then (iv) the flux $\mathbf{q}_i = B_s \nabla (\mathbf{\Pi}_0 u^*)_i$ is constant, (v) the source term ${q^*}$ is given by \eqref{eq-third}, and (vi) the model \eqref{eq-ell-modelf-global-pde} reads %% \begin{eqnarray} - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} \theta_s B_s \nabla (\mathbf{\Pi}_0 \nabla u^* )_i - \nabla \cdot ( \mathbf{{B}^*} \nabla u^* ) &=& {f^*}, \;\; \mathbf{x} \in \Omega. \label{eq-ell-modelf-1} \end{eqnarray} %% In both cases the solution $u^*$ is fully or partially decoupled from the local solution ${u_s^*}_i$ in the sense that the global problem can be solved independently of the local problem \eqref{eq-ell-modelf-block}--\eqref{eq-ell-modelf-bc}, which, in turn, can be solved given the boundary data. \end{corollary} %% We note that a continuous version of \eqref{eq-ell-modelf} with $\mathbf{\Pi}_0$ was considered in \mpcite{Horn97}{pp.145} where it was not noticed that $\mathbf{q}_i$ are null and ${q^*}=0$. The affine case was considered in \cite{spagnuolo01} but only the zero'th moment terms were computed and it was noted that $\mathbf{q}_i$ are constant but ${q^*}$ was not used. These facts are fundamental for the modeling pursued in this paper since the fluxes $\mathbf{q}_i$ are the prototypes of advective velocities used in the model of coupled flow-advection-diffusion. If $\mathbf{q}_i$ are zero, then no advection effects can be accounted for. If \eqref{eq-upscale-parabolic-obstacle} with ${\phi}^* = 0$ is used instead of \eqref{eq-ell-modelf-1}, then the effects of flow in inclusions is underpredicted. %% \begin{remark} \rm The global equation \eqref{eq-ell-modelf-1} should be compared to \eqref{eq-upscale-parabolic} with $\tilde{\phi} = 0$ as, in the continuous limit, it reads %% \begin{eqnarray*} - \nabla \cdot \left( (\theta_s \mathbf{B_{s}} + \mathbf{{B}^*} ) \nabla u^* \right) &=& f, \;\; \mathbf{x} \in \Omega. \label{eq-ell-upscale1} \end{eqnarray*} %% The effective coefficient $\theta_s \mathbf{B_{s}} + \mathbf{{B}^*} $ is close to but not the same as $\mathbf{\tilde{B}}$, however. The consequences for the flow problem will be discussed in Section~\ref{sec-newmodel}. \end{remark} %%%%%%%%%%%% \section{Discrete upscaled model with memory terms} \label{sec-convolution} In this section we show how to interpret the memory term ${q^*}=\Pi' \mathbf{q}$ arising in \eqref{eq-par-modelf} and how to partially decouple the system \eqref{eq-par-modelf}. The plan is to represent the solution ${u_s^*}_i$ to the cell problem \eqref{eq-par-modelf-block} as a linear functional acting on the boundary data given by $(\mathbf{\Pi} u^*)_i$ in \eqref{eq-par-modelf-bc}. With ${u_s^*}_i$ calculated, one computes its fluxes $q_i$; finally the values of $\Pi' \mathbf{q}$ to be inserted back to \eqref{eq-par-modelf-global-pde} follow. % These steps \[ u^* \mapsto \mathbf{\Pi} u^* \mapsto {u_s^*}_i \mapsto q_i \mapsto \Pi' \mathbf{q} \equiv {q^*} \] are a composition of linear functionals. Some are simple projections and extensions of polynomials local to $\Omega_{is}$ to functions on $\Omega$; these were defined in Section~\ref{sec-discrete}. The other functionals are Dirichlet-to-Neumann maps; these can be written using fundamental solutions or equivalently, as the time convolutions with auxiliary kernels depending only on the geometry of $\Omega_{is}$ and coefficients of the problem \eqref{eq-par-modelf-block}. Here we focus on these Dirichlet-to-Neumann calculations represented schematically by $\mathbf{\Pi} u^* \mapsto {u_s^*}_i \mapsto q_i $. The calculations depend on the choice of operator $\mathbf{\Pi}$: we first calculate $q_i$ when $\mathbf{\Pi}_0$ is used; next, we use $\mathbf{\Pi}_1$. The former model is related to the standard double-porosity model from \cite{ADH,Arb89}; the latter is related to the secondary diffusion model considered in \cite{ClarkShow94,CookShow95}. The calculations are done on a generic cell $\Omega_{0s}$, but it is easy to extend it to any $\Omega_{is}$. We focus on a generic self-adjoint parabolic model; calculations and definitions for the advection-diffusion model will be presented in Section~\ref{sec-newmodel}. Recall the notation and properties of the convolution product $\kappa \ast \lambda$ of any two functions $\kappa, \lambda \in L^1(0,T)$ defined as $ (\kappa \ast \lambda )(t) \equiv (\kappa(\cdot) \ast \lambda(\cdot))(t) \equiv \int_0^t \kappa(\tau) \lambda(t-\tau) d\tau, \;\; t \geq 0. $ %% The basic properties of this product include symmetry $\kappa \ast \lambda = \lambda \ast \kappa$ and the following differentiation %% $ \ddt{}(\kappa \ast \lambda) = \ddt{\kappa} \ast \lambda + \kappa(0)\lambda(t). = {\kappa} \ast \ddt{\lambda} + \kappa(t)\lambda(0). $ % %% Appropriate extensions are easily defined for vector valued functions $\kappa,\lambda \in L^1(0,T;X)$ where $X$ is some normed vector space, and similar properties hold, with $\ddt{}$ replaced by $\dtt{}$. %%%%%%%%%%%% \subsection{Memory terms from piecewise constant boundary conditions} \label{sec-convolution-primary} Here we provide the Dirichlet-to-Neumann calculations corresponding to $\mathbf{\Pi} _0u^* \mapsto {u_s^*}_0 \mapsto q_0 $, where ${u_s^*}_0$ solves \eqref{eq-par-modelf-block} with \eqref{eq-par-modelf-bc} and \eqref{eq-par-modelf-inits} on a generic cell $\Omega_{0s}$. Similar calculations were done in \cite{Arb89,HornSho90,P92}. Consider a representative solution $r^0=r^0(\mathbf{y},t)$ which solves % \begin{eqnarray} \label{eq-def-r0} \left\{ \begin{array}{lll} \phi_s \frac{\partial r^0}{\partial t} - \nabla \cdot ( \mathbf{B_{s}} \nabla r^0) &= 0, \;\; &\mathbf{y} \in \Omega_{0s} \\ r^0(\mathbf{y},0) &= 0, \;\; &\mathbf{y} \in \Omega_{0s} \\ r^0(\mathbf{y},t) &= 1, \;\; &\mathbf{y} \in \Gamma_0 \end{array} \right. \end{eqnarray} %%% It is also convenient to consider $r(\mathbf{y},t)=1-r^0(\mathbf{y},t)$ which solves %% \begin{eqnarray} \label{eq-def-r} \left\{ \begin{array}{lll} \phi_s \frac{\partial r}{\partial t} - \nabla \cdot ( \mathbf{B_{s}} \nabla r) &= 0, \;\; &\mathbf{y} \in \Omega_{0s} \\ r(\mathbf{y},0) &= 1, \;\; &\mathbf{y} \in \Omega_{0s} \\ r(\mathbf{y},t) &= 0, \;\; &\mathbf{y} \in \Gamma_0. \end{array} \right. \end{eqnarray} %% Next we define %% \begin{eqnarray} \label{eq-kernel00} \mathcal{T}^{00} (t) \equiv \frac{1}{\abs{\Omega_0}} \int_{\Omega_{0s}} \phi_s\frac{\partial r^0(\mathbf{y},t)}{\partial t} dA \end{eqnarray} %% This is the first of the kernels to be used in the sequel. Note that % \begin{eqnarray} \label{eq-kernel00-r} \mathcal{T}^{00} (t) = - \frac{1}{\abs{\Omega_0}} \int_{\Omega_{0s}} \phi_s\frac{\partial r(\mathbf{y},t)}{\partial t} dA. \end{eqnarray} %%%%%%%%%%%%% \begin{proposition} \label{prop-conv-0} Let $\mathbf{\Pi}=\mathbf{\Pi}_0$ in \eqref{eq-par-modelf} and let there be an initial equilibrium, that is, let %% \begin{eqnarray} \label{eq-equilibrium} u_i^0 = (\mathbf{\Pi}_0 u^0)_i = u^0 \equiv const, \; \forall i \end{eqnarray} %% in \eqref{eq-par-modelf-initf} and \eqref{eq-par-modelf-inits}. Let also the assumption on periodic geometry hold in the sense that all inclusions have congruent geometry and equal coefficients: %% \begin{eqnarray} \label{eq-congr} \Omega_{0s} \cong \Omega_{is}, \;\;{\rm\ and\ } \phi_{is}=\phi_s, \mathbf{B}_{is}=B_s \; \forall i=1, \dots,N_{\rm incl} \end{eqnarray} %% Then the equation \eqref{eq-par-modelf-global-pde} can be written in the convolution form as follows %%% \begin{eqnarray} {\phi}^* \frac{\partial u^*}{\partial t} + \sum_i \hat{\chi}_i(\mathbf{x}) \mathcal{T}^{00}(t) \ast \ddt{(\mathbf{\Pi}_0 u^*)_i} - \nabla \cdot ( \mathbf{{B}^*} \nabla u^* ) &=& 0, \;\; \mathbf{x} \in \Omega. \label{eq-par-convolution} \end{eqnarray} %% \end{proposition} \begin{proof} This is obtained by linearity from the following calculations. Let $u_0^0 \equiv const$ be given and $A^0:[0,\infty)\rightarrow \mathbb{R}$ be a given differentiable function, continuous at $0$. Define %% \begin{eqnarray} {u_s^*}_0(\mathbf{y},t) \equiv \ddt{A^0(t)} \ast r^0(\mathbf{y},t) +A^0(0)r^0(\mathbf{y},t) + u_0^0 (1- r^0(\mathbf{y},t)). \label{eq-convolution-u-0} \end{eqnarray} %% We note in passing that ${u_s^*}_0(\mathbf{y},t)= A^0(\cdot)\ast \dtt{r^0(\mathbf{y},\cdot)} + u_0^0 r(\mathbf{y},t)$, which follows by differentiating convolutions from $A^0 \ast \dtt{r^0} = A^0 \ast \dtt{r^0} +A^0(t)r^0(\mathbf{y},0) = \ddt{} (A^0 \ast r^0) = \ddt{} (r^0 \ast A^0) = \ddt{A^0} \ast r^0 + A^0(0) r^0(t) $. We easily verify that ${u_s^*}_0$ satisfies % \begin{subequations} \label{eq-block-u} \begin{eqnarray} \phi_s \frac{\partial {u_s^*}_0}{\partial t} - \nabla \cdot ( \mathbf{B_{s}} \nabla {u_s^*}_0) &=& 0, \;\; \mathbf{y} \in \Omega_{0s} \label{eq-block-u-pde} \\ {u_s^*}_0(\mathbf{y},0) \equiv const &=& u_0^0, \;\; \mathbf{y} \in \Omega_{0s} \label{eq-block-u-init} \\ {u_s^*}_0(\mathbf{y},t) &=& A^0(t), \;\; \mathbf{y} \in \Gamma_0 \label{eq-block-u-bc} \end{eqnarray} \end{subequations} % Indeed, we calculate the first term in \eqref{eq-block-u-pde} %% \begin{multline*} \phi_s \dtt{{u_s^*}_0(\mathbf{y},t)} = \phi_s \left( \ddt{A_0^0}\ast \dtt{r^0(\mathbf{y},t) } + \ddt{A_0^0} (t) r^0(\mathbf{y},0) + (A_0^0(0)-u_0^0) \dtt{r^0(\mathbf{y},t)} \right) \\ = \ddt{A_0^0}\ast \phi_s \dtt{r^0(\mathbf{y},t) } + (A_0^0(0)-u_0^0)\phi_s \dtt{r^0(\mathbf{y},t)}, \end{multline*} %% as well as the second term %% \[ - \nabla \cdot B_s \nabla {u_s^*}_0(\mathbf{y},t) = \ddt{A^0} \ast \left(-\nabla \cdot B_s \nabla r^0(\mathbf{y},t)\right) -(A_0^0(0) -u_0^0)\left(\nabla \cdot B_s \nabla r^0(\mathbf{y},t)\right). \] % Combine these and use \eqref{eq-def-r0} to see that \eqref{eq-block-u-pde} is satisfied. The initial condition \eqref{eq-block-u-init} is trivially satisfied when $t=0$ is used in \eqref{eq-convolution-u-0}. Finally, the boundary condition \eqref{eq-block-u-bc} follows from \eqref{eq-def-r0}, by %% \begin{multline*} {u_s^*}_0(\mathbf{y},t) \vert_{\Gamma_0} = \left( \ddt{A_0^0(t)} \ast r^0(\mathbf{y},t) +A_0^0(0)r^0(\mathbf{y},t) + u_0^0(1-r^0(\mathbf{y},t))\right) \vert_{\Gamma_0} \\ = \ddt{A_0^0(t)} \ast 1 +A_0^0(0) = A_0^0(t)-A_0^0(0) + A_0^0(0)=A_0^0(t). \end{multline*} % Next, compute the total flux out of $\Omega_{0s}$ by the divergence theorem, %% \begin{multline*} \int_{\Gamma_0} B_s \nabla {u_s^*}_0 \cdot \mathbf{n} = \int_{\Omega_{0s}} \phi_s \dtt{{u_s^*}_0} dA \\ = \ddt{A_0^0} \ast \left(\int_{\Omega_{0s}} {\phi_s \frac{\partial r^0(\mathbf{y},t)}{\partial t} } dA\right) +\left(A_0^0(0)-u_0^0\right) \int_{\Omega_{0s}} \phi_s \frac{\partial r^0(\mathbf{y},t)}{\partial t} dA \\ = \ddt{A_0^0} \ast \left( - \int_{\Omega_{0s}} {\phi_s \frac{\partial r(\mathbf{y},t)}{\partial t} } dA\right) %\\ +\left(A_0^0(0)-u_0^0\right) \left( - \int_{\Omega_{0s}} \phi_s \frac{\partial r(\mathbf{y},t)}{\partial t} dA\right). \end{multline*} %% Now use \eqref{eq-kernel00} to conclude that $ \abs{\Omega_0} \mathbf{M}^0 (B_s \nabla {u_s^*}_0 \cdot \mathbf{n}) = \int_{\Gamma_0} B_s \nabla {u_s^*}_0 \cdot\mathbf{n} $ and %% \begin{eqnarray} \label{eq-conv-flux} \mathbf{M}^0 (B_s \nabla {u_s^*}_0 \cdot \mathbf{n}) = \ddt{A^0(\cdot)} \ast \mathcal{T}^{00} (\cdot) + \left(A_0^0(0)-u_0^0\right) \mathcal{T}^{00}(t). \end{eqnarray} %% We have thus computed the flux compatible with the solution to \eqref{eq-block-u}. Now, for each $i$, we apply analogous calculations to the solutions to \eqref{eq-par-model-block} with boundary condition \eqref{eq-par-model-bc} and initial condition \eqref{eq-par-model-inits} over $\Omega_{is}$. The boundary condition analogous to the one in \eqref{eq-block-u-bc} is $A_i^0(t) = (\mathbf{\Pi}_0 u^*)_i$, and the initial condition analogous to \eqref{eq-block-u-init} is given by $u_i^0(t)$. % We get $M^0_i(\mathbf{q})$ with $q_i = B_s \nabla {u_s^*}_i \cdot \mathbf{n}$ as in \eqref{eq-conv-flux}, %% \[ M_i^0 (\mathbf{B_{s}} \nabla {u_s^*}_i(\mathbf{s},t) \cdot \mathbf{n}) = \mathcal{T}^{00}(t) \ast \ddt{(\mathbf{\Pi}_0 u^*)_i} + (A^0_i(0)-u_i^0) \mathcal{T}^{00}(t) . \] Finally, by \eqref{eq-equilibrium} we have initial equilibrium so that the last term above vanishes. Substituting in \eqref{eq-par-model-global-pde} yields %% \[ {q^*}(\mathbf{x},t) = \Pi' \mathbf{q} = \sum_i \hat{\chi}_i(\mathbf{x}) M_i^0 (\mathbf{B_{s}} \nabla {u_s^*}_i(\mathbf{s},t) \cdot \mathbf{n}) = \sum_i \hat{\chi}_i(\mathbf{x}) \mathcal{T}^{00}(t) \ast \ddt{(\mathbf{\Pi}_0 u^*)_i}. \] % \end{proof} \begin{remark} \rm The character of $\mathcal{T}^{00}(t)$ is easily understood from the characterization \eqref{eq-kernel00-r}: $\mathcal{T}^{00}$ is monotone decreasing, unbounded at $0$, and positive in the sense pursued in \cite{McleanThomee}. Its singularity at $0$ is weak in the sense that the (improper) integral $\int_0^t \mathcal{T}^{00}(s) ds$ is finite. One can see that for small times $t$, the kernel $\mathcal{T}^{00}(t)$ behaves like $t^{-\alpha}$ with $\alpha = 1/2$ \cite{CarslawJaeger}; details are provided in \cite{P92}. \end{remark} \begin{remark} \rm The assumptions \eqref{eq-equilibrium} and \eqref{eq-congr} can be easily relaxed; the resulting global equation has additional terms corresponding to additional tailing effects. In addition, for each $i=1, \ldots N_{\rm incl}$, a different kernel $\mathcal{T}^{00}_i$ is constructed. \end{remark} The system \eqref{eq-par-convolution} is a single integro-differential equation which is formally equivalent to the coupled system \eqref{eq-par-modelf}. Theoretically, the kernel $\mathcal{T}^{00}$ can be pre-computed analytically for simple geometries \cite{CarslawJaeger,P92} or numerically for more general cases \cite{P92}. This direction was also pursued in \cite{Jaffre02}. From an analytical and modeling point of view the single equation \eqref{eq-par-convolution} is very attractive. However, in a computational realization, the coupled form \eqref{eq-par-modelf} may be preferred due to the following issues. The (mild) difficulties in direct discretization of \eqref{eq-par-convolution} arise due to singularity of $\mathcal{T}^{00}(t)$ at $0$. Somewhat more limiting are the long-term memory effects which require storing all history of $u^*$ if \eqref{eq-par-convolution} is used. The latter can be alleviated if the history is truncated, as discussed theoretically in \cite{MclTho,P95c} and as is frequently done in applications \cite{Hagg95,Hagg00,Hagg01}. %%%%%%%%%% \subsection{Memory terms from piecewise affine boundary conditions} \label{sec-convolution-secondary} Consider now the model \eqref{eq-par-modelf} in which we use $\mathbf{\Pi}_1,\Pi_1'$. We provide the calculations of the Dirichlet-Neumann map $\mathbf{\Pi} _1u^* \mapsto {u_s^*}_0 \mapsto q_0 $, and derive the effective model equivalent to \eqref{eq-par-modelf} in a convolution form; this ``secondary diffusion'' version was considered in \cite{CookShow95}. As in Section~\ref{sec-convolution-primary}, we consider solutions to the cell problem % \begin{subequations} \label{eq-block-u-2} \begin{eqnarray} \phi_s \frac{\partial {u_s^*}_0}{\partial t} - \nabla \cdot ( \mathbf{B_{s}} \nabla {u_s^*}_0) &=& 0, \label{eq-block-u-pde-2} \\ {u_s^*}_0(\mathbf{y},0) &=& u_0^0, \;\; \mathbf{y} \in \Omega_{0s} , \label{eq-block-u-init-2} \\ {u_s^*}_0(\mathbf{y},t) = A_0^0(t) &+& \sum_{k=1}^2 A_0^k(t) (\mathbf{y}-\mathbf{x}_0^c)_k, \;\; \mathbf{y} \in \Gamma_0\,, \label{eq-block-u-bc-2} \end{eqnarray} % \end{subequations} % subject to an affine \eqref{eq-block-u-bc-2} rather than a constant boundary condition as in \eqref{eq-block-u-bc}. The coefficients $A_0^1(t), \; A_0^2(t)$ distinguish the model \eqref{eq-block-u-2} from \eqref{eq-block-u}. For these we define additional auxiliary functions, for $k=1,2$ % \begin{eqnarray} \label{eq-def-rk} \left\{ \begin{array}{lll} \phi_s \frac{\partial r^k}{\partial t} - \nabla \cdot ( \mathbf{B_{s}} \nabla r^k) &= 0, \;\; &\mathbf{y} \in \Omega_{0s} , \\ r^k(\mathbf{y},0) &= 0, \;\; &\mathbf{y} \in \Omega_{0s} , \\ r^k(\mathbf{y},t) &= (\mathbf{y}-\mathbf{x}^c_0)_k, \;\; &\mathbf{y} \in \Gamma_0. \end{array} \right. \end{eqnarray} %% %% \begin{lemma} %% Let $A_0^1(0)=A_0^2(0)=0$ hold and let $A_0^0=u_0^0$. Define %% \begin{multline} {u_s^*}_0(\mathbf{y},t) = \ddt{A^0(\cdot)} \ast {r^0(\mathbf{y},\cdot)} + A_0^0(0) r^0(\mathbf{y},t) + u_0^0 (1-r^0(\mathbf{y},t)) \\ + \sum_{k=1}^2 A^k(\cdot) \ast \frac{\partial r^k}{\partial t} (\mathbf{y},\cdot ) \,. \; \mathbf{y} \in \Omega_0. \label{eq-convolution-u-1} \end{multline} Then ${u_s^*}_0$ solves \eqref{eq-block-u-2}. \label{lem-represent-u} \end{lemma} %% \begin{proof} The proof follows from calculations similar to those in Proposition \eqref{prop-conv-0}. We verify the additional terms coming from the last part of \eqref{eq-convolution-u-1}. In fact, by $A_0^1(0)=A_0^2(0)=0$ and $r^k(\mathbf{y},0)\equiv 0$ we have $\ddt{} (A^k \ast r^k) = \ddt{A^k} \ast r^k = A^k \ast \dtt{r^k}$ for $k=1,2$. We compute, for $\mathbf{y} \in \Omega_{0s}$, using $A^0(0)=u^0_0$ %% \begin{eqnarray*} \dtt{{u_s^*}_0}(\mathbf{y},t) &=& \sum_{k=0}^2 \ddt{A^k}(\cdot) \ast \frac{\partial r^k}{\partial t} (\mathbf{y},\cdot ), \label{eq-convolution-u-1-dt} \\ \mathbf{B_{s}} \nabla {u_s^*}_0(\mathbf{y},t) &=& \sum_{k=0}^2 A^k(\cdot) \ast \mathbf{B_{s}} \nabla \frac{\partial r^k (\mathbf{y},\cdot )}{\partial t} = \sum_{k=0}^2 \ddt{A^k} (\cdot) \ast \mathbf{B_{s}} \nabla r^k (\mathbf{y},\cdot ) \\&=& A^0 \ast \dtt{} \mathbf{B_{s}} \nabla r^0 - u^0_0 \mathbf{B_{s}} \nabla r^0 + \sum_{k=1}^2 A^k(\cdot) \ast \mathbf{B_{s}} \nabla \frac{\partial r^k (\mathbf{y},\cdot )}{\partial t}, \label{eq-convolution-u-nabla} \end{eqnarray*} %% and immediately verify that $u^*$ defined by \eqref{eq-convolution-u-1} satisfies the PDE, the boundary and initial conditions of \eqref{eq-block-u-2}. \end{proof} Now we follow similar steps as in Section~\ref{sec-convolution-primary} to represent ${q^*}=\Pi_1' \mathbf{q}$. We use kernels arising from various averages of $r^k$. First, we use the averages of rate of change in time %% \begin{eqnarray} \label{eq-kernelk0} \mathcal{T}^{k0}(t) \equiv \frac{1}{\abs{\Omega_0}} \int_{\Omega_{0s}} \phi_{s} \frac{\partial r^k}{\partial t} (\mathbf{y},t)\,dA, \quad k=0,1,2 \end{eqnarray} %% where $\mathcal{T}^{00}$ defined previously in \eqref{eq-kernel00} is included for completeness. Next, the kernels $\mathcal{T}^{k1},\mathcal{T}^{k2}$ arising from the first moments of $r^k,k=0,1,2$ are defined by %% \begin{eqnarray} \label{eq-kernel0k} \mathcal{T}^{kj}(t) \equiv \frac{1}{\abs{\Omega_0}} \int_{\Omega_{0s}} \phi_s \frac{\partial r^k}{\partial t} (\mathbf{y},t)(\mathbf{y}- (\mathbf{x}^C_0))_j\,dA, \quad j=1,2;\; k=0,1,2. \end{eqnarray} %% Finally, for each $r^k,k=0,1,2$ we set %% \begin{eqnarray} \label{eq-kernelkk} \mathbf{S}^k(t) \equiv (S^{k1},S^{k2}) \equiv \left[ \begin{array}{l}S^{k1}\\S^{k2}\end{array} \right] \equiv \frac{1}{\abs{\Omega_0}} \int_{\Omega_{0s}} \mathbf{B_{s}} \nabla r^k (\mathbf{y},t)\,dA. \end{eqnarray} %% In summary, we have defined the total of fifteen geometry-based and time-dependent kernels: nine zero'th and first order moments $\mathcal{T}^{k0},\mathcal{T}^{k1},\mathcal{T}^{k2}$ of $r^k,k=0,1,2$ and six averages $S^{k1},S^{k2}$ for $k=0,1,2$ of their gradients $\nabla r^k$. These are used to express ${q^*}$ in the upscaled model \eqref{eq-par-model}. If \eqref{eq-congr} does not hold, then these kernels can be defined separately for each $i$. We note that many of these kernels may vanish due to symmetry; see Remark~\ref{rem-symmetry}. We can now represent the solution ${u_s^*}_i$, for each $i$, to \eqref{eq-par-model-block} with \eqref{eq-par-model-bc} and \eqref{eq-par-model-inits}, as it was done in Lemma~\ref{lem-represent-u} for the solution ${u_s^*}_0$ to \eqref{eq-block-u-2}. We use boundary conditions expressed by $A_i^k,k=0,1,2$. Clearly $A_i^k,k=0,1,2$ and thus ${u_s^*}_i$ vary with $i$. %% Next we compute ${q^*}=\Pi_1' \mathbf{q}$ where $\mathbf{q}=(q_i)_{i=1}^{N_{\rm incl}}, \;q_i=\mathbf{q}_i \cdot \mathbf{n}_i = \mathbf{B_{s}} \nabla {u_s^*}_i \cdot \mathbf{n}_i$. % By Lemma~\ref{lem-ppi1}, the following terms arise as in \eqref{eq-moment1div}: % \[\mathbf{\Pi}_1' \mathbf{q}= \sum_i \hat{\chi}_i(\mathbf{x}) I_i - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} II_i - \nabla \cdot \sum_i {\hat{\chi}_i(\mathbf{x})} III_i\,. \] % Here $I$ is a scalar and $II,III \in \mathbb{R}^2$ are vectors; we write $II=(II^1,II^2)$ and $III=(III^1,III^2)$. These terms are as follows % \begin{eqnarray*} I_i & \equiv& \frac{1}{\abs{\Omega_0}} \int_{\Omega_{is}}( \nabla \cdot \mathbf{B_{s}} \nabla {u_s^*}_i) dA, \\ II_i & \equiv& \frac{1}{\abs{\Omega_0}} \int_{\Omega_{is}} (\nabla \cdot \mathbf{B_{s}} \nabla {u_s^*}_i ) (\mathbf{y} - {\mathbf{x}}^C_i) dA, \\ III_i & \equiv& \frac{1}{\abs{\Omega_{0s}}} \int_{\Omega_{is}} (\mathbf{B_{s}} \nabla {u_s^*}_i ) dA. \end{eqnarray*} %% The zero'th moments of rate of change of mass content are % \begin{eqnarray*} I_i = \frac{1}{\abs{\Omega_0}} \int_{\Omega_{0s}} \phi_s \dtt{{u_s^*}_i}(\mathbf{y},t)\,dA = \sum_{k=0}^2 \mathcal{T}_i^{k0} \ast \ddt{A_i^k}(\cdot). \end{eqnarray*} %% Next we calculate the { first moment} of mass content in $\Omega_{is}$ %% \begin{multline*} II_i = \frac{1}{\abs{\Omega_0}} \int_{\Omega_{is}} \phi_{s} \dtt{{u_s^*}_i}(\mathbf{y},t)(\mathbf{y} - {\mathbf{x}}^C_i)\,dA \\ = \sum_{k=0}^2 \left( \mathcal{T}_i^{k1}(\cdot) \ast \ddt{A_i^k}(\cdot),\mathcal{T}_i^{k2}(\cdot) \ast \ddt{A_i^k}(\cdot)\right) . \end{multline*} %% In the last step we handle $III_i$: %% \begin{eqnarray*} III_i = \frac{1}{\abs{\Omega_0}} \int_{\Omega_{is}} (\mathbf{B_{s}} \nabla {u_s^*}_i(\mathbf{y},t) )\,dA = \sum_{k=0}^2 \mathbf{S}_i^k(\cdot) \ast \ddt{A_i^k}(\cdot). \end{eqnarray*} %% Now we substitute explicitly $A_i^k, k=0,1,2$ in \eqref{eq-block-u-bc-2} via $\mathbf{\Pi}_1 u^*$ as in \eqref{eq-par-model-bc} %% \begin{eqnarray} \label{eq-a0} A_i^0(t) &=& (\mathbf{\Pi}_1 u^*)_i = \frac{1}{\snorm{\Omega_i}{}} \int_{\Omega_{is}} u^*(\mathbf{y},t)\,dA\,, \\ \label{eq-ak} A_i^k(t) &=& (\mathbf{\Pi}_1 \partial_k u^*)_i =\frac{1}{\snorm{\Omega_i}{}} \int_{\Omega_{is}} \partial_k u^*(\mathbf{y},t)\,dA\,, k=1,2 \end{eqnarray} %% and complete the calculation %% \begin{multline} \label{eq-represent-q} \Pi_1' \mathbf{q} = \sum_{i} \hat{\chi}_i(\mathbf{x}) \sum_{k=0}^2 \mathcal{T}_i^{k0} \ast \ddt{A_i^k}(\cdot) \\ - \nabla \cdot \sum_{i} \hat{\chi}_i(\mathbf{x}) \sum_{k=0}^2 (\mathcal{T}_i^{k1}(\cdot),\mathcal{T}_i^{k2}(\cdot)) \ast \ddt{A_i^k}(\cdot) \\ - \nabla \cdot \left( \sum_{i} \hat{\chi}_i(\mathbf{x}) \sum_{k=0}^2 \mathbf{S}_i^k(\cdot)\ast \ddt{A_i^k}(\cdot)) \right) \,. \end{multline} %% The final step of representing $\Pi_1' \mathbf{q}$ follows after we take advantage of \eqref{eq-a0},\eqref{eq-ak}, and insert \eqref{eq-represent-q} to \eqref{eq-par-model}. This completes the proof of the next Proposition. \begin{proposition} \label{prop-pi1} Let the assumptions of Lemma~\ref{lem-represent-u} hold. Let also \eqref{eq-congr} hold so we can supress the index $i$ on each of the kernels. Then the global PDE \eqref{eq-par-modelf} using \eqref{eq-represent-q} takes the convolution form %% \begin{multline} \label{eq-gl-conv-pde} % {\phi}^* \dtt{u^*(\mathbf{x},t)} % + \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} \mathcal{T}^{00}(\cdot) \ast \ddt{(\mathbf{\Pi}_0 u^*(,\cdot))_i} + \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} \sum_{k=1}^2 \mathcal{T}^{k0}(\cdot) \ast \ddt{ (\mathbf{\Pi}_0 \partial_ku^*(\cdot))_i} \\ - \nabla \cdot \left( % \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} [ \mathcal{T}^{01}_i(\cdot ),\mathcal{T}^{02}_i(\cdot)]^T \ast \ddt{ (\mathbf{\Pi}_0 u^*(\cdot))_i}\right. \\ +\left. \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} \sum_{k=1}^2 [ \mathcal{T}^{k1}_i(\cdot ),\mathcal{T}^{k2}_i(\cdot)]^T \ast \ddt{ (\mathbf{\Pi}_0 \partial_k u^*(\cdot))_i}\right. \\ \left. + \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} [S^{01}_i(\cdot),S^{02}_i(\cdot)]^T \ast \ddt{ (\mathbf{\Pi}_0 u^*(\cdot))_i } \right. \\ + \left. \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} \sum_{k=0}^2 [S^{k1}_i(\cdot),S^{k2}_i(\cdot)]^T \ast \ddt{ (\mathbf{\Pi}_0\partial_k u^*(\cdot))_i } \right) % \\ - \nabla \cdot \left( \mathbf{{B}^*} \nabla u^*(\mathbf{x},t) \right) = 0, \mathbf{x} \in \Omega,\ t > 0. \end{multline} %% \end{proposition} \begin{remark} \label{rem-nonpicky} \rm Let $N_{\rm incl}$ be large. Then, by \eqref{eq-limit} we can ``formally'' let \newline $\hat{\chi}_i(\mathbf{x}) (\mathbf{\Pi}_0 \partial_k u^*)_i \to \partial_k u^*(\mathbf{x},t), k=0,1,2$. Also, the weak derivative can be interpreted as %\newline $\nabla (\sum_i \hat{\chi}_i(\mathbf{x}) (\mathbf{\Pi}_0 u^*)_i) \to \nabla u^*(\mathbf{x},t)$. % With these informal limits, the structure of the limiting model is % \begin{multline} {\phi}^* \dtt{u^*} + \mathcal{T}^{00} \ast \dtt{u^*} \\ + (\mathcal{T}^{10},\mathcal{T}^{20}) \ast \nabla \dtt{u^*} - \nabla \cdot \left( (\mathcal{T}^{01},\mathcal{T}^{02} ) \ast \dtt{u^*} \right) - \nabla \cdot \left( (S^{01},S^{02}) \ast \dtt{u^*} \right) \\% -\nabla \cdot \left( \left[ \begin{array}{ll} \mathcal{T}^{11}&\mathcal{T}^{12}\\ \mathcal{T}^{21}&\mathcal{T}^{22} \end{array}\right] \ast \nabla \dtt{u^*} \right) %\\ - \nabla \cdot \left( \left[ \begin{array}{ll} S^{11}&S^{12}\\ S^{21}&S^{22} \end{array}\right] \ast \nabla \dtt{u^*} \right) \\ - \nabla \cdot \left( \mathbf{{B}^*} \nabla u^* \right) = 0 \end{multline} % or, after we collect like-terms % \begin{eqnarray} \label{eq-nonpicky} {\phi}^* \dtt{u^*} + \mathcal{T}^{00} \ast \dtt{u^*} + \mathbf{M} \ast \nabla \dtt{u^*} -\nabla \cdot \left( \mathcal{M} \ast \nabla \dtt{u^*} \right) - \nabla \cdot \left( \mathbf{{B}^*} \nabla u^* \right) = 0 \end{eqnarray} % where $\mathbf{M}:(0,\infty) \rightarrow \mathbb{R}^2, \mathcal{M}:(0,\infty) \rightarrow \mathbb{R}^{2 \times 2}$ are time dependent vector and matrix valued memory kernels. \end{remark} \begin{remark} \label{rem-symmetry} \rm Let $\Omega_{is}$ be circular or square inclusions. Then by symmetry, the kernels $ \mathcal{T}^{10}, \mathcal{T}^{20},$ $ \mathcal{T}^{01}, \mathcal{T}^{02}, \mathcal{T}^{12},\mathcal{T}^{21}, S^0$ vanish. In addition, $ \left[ \begin{array}{ll} \mathcal{T}^{11}&\mathcal{T}^{12}\\ \mathcal{T}^{21}&\mathcal{T}^{22} \end{array}\right]$ and $[\mathbf{S}^1,\mathbf{S}^2] $ are diagonal. Then the model \eqref{eq-nonpicky} becomes the secondary diffusion model from \cite{ClarkShow94} % \begin{eqnarray} \label{eq-nonpicky2} {\phi}^* \dtt{u^*} + \mathcal{T}^{00} \ast \dtt{u^*} -\nabla \cdot \left( \mathcal{M} \ast \nabla \dtt{u^*} \right) %\\ - \nabla \cdot \left( \mathbf{{B}^*} \nabla u^* \right) = 0, \end{eqnarray} and $\mathcal{M}$ is diagonal. \end{remark} In many physical situations the symmetries mentioned in Remark~\ref{rem-symmetry} hold and thereby the terms associated with $\mathcal{M}$ in the effective model~\eqref{eq-nonpicky2} are most significant at most time scales. Numerical evidence assessing practical importance of secondary diffusion terms versus all other terms will be presented elsewhere. We stress that such a model arises via upscaling of a self-adjoint parabolic problem. On the other hand, in non-self adjoint problems, the terms associated with off-diagonal kernels will not vanish. For example, problems with first order terms such as those in advection-diffusion-dispersion problems to be discussed in Section~\ref{sec-newmodel}, will not simplify to \eqref{eq-nonpicky2}. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Upscaling the coupled flow-advection-diffusion model} \label{sec-newmodel} Now we return to the {\em flow-advection-diffusion system} \eqref{eq-flow}--\eqref{eq-transport}, or \eqref{eq-adv-dif-flo-int}--\eqref{eq-adv-dif-tra-int}. % The results in Sections~\ref{sec-review}, \ref{sec-discrete}, \ref{sec-convolution}, do not apply directly to \eqref{eq-adv-dif-tra-int} due to a) non-symmetry due to advection, and due to b) the coupled nature of the system: the scales of diffusion, advection, and dispersion in \eqref{eq-adv-dif-tra-int} are coupled to the scales of flow in \eqref{eq-adv-dif-flo-int}. We first recall the results without any scaling as well as in the ``obstacle'' limit for \eqref{eq-adv-dif-flo-int}--\eqref{eq-adv-dif-tra-int}. Next we discuss the consequences of $\epsilon^2$-scaling to \eqref{eq-adv-dif-flo-int}--\eqref{eq-adv-dif-tra-int}, that is, classical double porosity approaches which correspond to using $\mathbf{\Pi}_0$, but in which the advection and dispersion effects are lost. Then we propose the affine approximations associated with operator $\mathbf{\Pi}_1$ instead. The final step is to represent the various memory terms arising via a Dirichlet-to-Neumann map in convolution form as in Section~\ref{sec-convolution}. %%%%%%%%%%%%%% \subsection{Effective coupled model by classical upscaling and obstacle limit} \label{sec-newmodel-classical} The solution of the exact system \eqref{eq-adv-dif-flo-int} is denoted by $p(\mathbf{x})$, $\mathbf{v}(\mathbf{x})$, $c(\mathbf{x},t)$. Using classical upscaling techniques these can be approximated by their local averages: the respective upscaled functions $\tilde{p}$, $\mathbf{\tilde{v}}$, and $\tilde{c}$ satisfy the approximating upscaled model % \begin{subequations} \label{eq-adv-dif-ups} \begin{eqnarray} \label{eq-divfree} \nabla \cdot \mathbf{\tilde{v}} &=& 0, \;\; \mathbf{x} \in \Omega, \\ \mathbf{\tilde{v}} &=& - \mathbf{\tilde{K}} \nabla \tilde{p}, \end{eqnarray} %% in which the effective conductivity $\mathbf{\tilde{K}} = \mathbf{\tilde{K}}(K_f,K_s)$ is given analogously to $\mathbf{\tilde{B}}$ as in \eqref{eq-upscaleBB}. The effective $\mathbf{\tilde{v}}$ is divergence-free by \eqref{eq-divfree} and we can use it directly in the upscaled transport model %% \begin{eqnarray} \tilde{\phi} \frac{\partial \tilde{c}}{\partial t} + \nabla \cdot ( \mathbf{\tilde{v}} \tilde{c} - \tilde{\mathbf{D}} \nabla \tilde{c}) = 0, \;\; \mathbf{x} \in \Omega, \end{eqnarray} \end{subequations} % where the effective {constant coefficients} $\tilde{\phi}=\tilde{\phi}(\phi_f,\phi_s)$ and $\tilde{\mathbf{D}}=\tilde{\mathbf{D}}(\mathbf{D_{f}},\mathbf{D_s})$ are computed with \eqref{eq-def-tphi} and \eqref{eq-upscaleB}, respectively. % The corresponding theoretical results can be found in \mpcite{Horn97}{pp 8-12, 243-246}, \mpcite{Allaire92}{Section 2,Thm 2.2}. Discussion of first order terms is in \mpcite{BLP}{pp 181-185} and \mpcite{JKO}{p 31}. As mentioned before in Section~\ref{sec-double-porosity}, this model is good only for the low contrast case and does not capture the tailing effects associated with storage in $\Omega_s$. Therefore we need to use the double porosity concept which is equivalent to the obstacle problem with a memory term. The obstacle case for the coupled system in which the blocks or inclusions are {\em impermeable} is obtained from the preceding case by setting formally $\mathbf{K}_s = 0,\phi_s = 0,\mathbf{D}_s = 0$, so we obtain the effective coefficients $\mathbf{\tilde{K}^0},\tilde{\phi}^0,\mathbf{\tilde{D}^0}$, and take note that $\Omega_f$ needs to be connected. Note that the upscaled unknowns $\mathbf{\tilde{v}}^0(x),\ \tilde{p}^0(x),\ \tilde{c}^0(x,t)$ are defined on all of $\Omega$. They are obtained as the solution of the system % \begin{subequations} \label{eq-adv-dif-obs} \begin{eqnarray} \nabla \cdot {\mathbf{\tilde{v}}^0} &=& 0, \;\; \mathbf{x} \in \Omega, \\ {\mathbf{\tilde{v}}^0} &=& - \mathbf{\tilde{K}^0} \nabla \tilde{p}^0, \;\; \mathbf{x} \in \Omega, \\ \tilde{\phi}^0 \frac{\partial \tilde{c}^0}{\partial t} + \nabla \cdot ( {\mathbf{\tilde{v}}^0} \tilde{c}^0 - \mathbf{\tilde{D}^0} \nabla \tilde{c}^0) &=& 0, \;\; \mathbf{x} \in \Omega, \end{eqnarray} \end{subequations} %% with appropriate boundary and initial conditions. Discussions of this case can be found in \mpcite{Horn97}{pp 13-16}, and \mpcite{Allaire92}{Thm 2.7}. Of course we hardly have $K_s=0$, this model is only auxiliary, and we are now ready to account for additional storage in $\Omega_s$. %%%%%%%%%%%%%%%%%% \subsection{Effective coupled model by discrete double porosity approach} \label{sec-newmodel-double-porosity} Now we follow ideas from Section~\ref{sec-double-porosity} and derive a general discrete upscaled model for \eqref{eq-adv-dif-flo-int}--\eqref{eq-adv-dif-tra-int}. First we revisit the discrete double-porosity model for the flow, following Section~\ref{sec-discrete-elliptic} and Corollary~\ref{cor-ell}. % We get the system for the flow %% \begin{subequations} \label{eq-par-flow-adv-dif} \begin{eqnarray} \nabla \cdot {\overline{\mathbf{v}^*}} \equiv - \Pi' \left( \ivec{{\mathbf{v}_s^*}_i \cdot \mathbf{n}_i} \right) + \nabla \cdot {\mathbf{v}^*} &=& 0, \\ {\mathbf{v}^*} &=& - \mathbf{{K}^*} \nabla {p^*}, \; \; \mathbf{x} \in \Omega , \label{eq-par-f-global-pde} %% \\ \nabla \cdot {\mathbf{v}_s^*}_i &=& 0, \; i=1,\ldots \in N_{\rm incl} \\ {\mathbf{v}_s^*}_i &=& - \mathbf{K}_s \nabla {p_s^*}_i , \;\; \mathbf{y} \in \Omega_{is} \label{eq-flow2-upscaled-loc} \\ % {p_s^*}_i \vert_{\Gamma_i} &=& (\mathbf{\Pi}({p^*}))_i\,, %% \label{eq-flow2-upscaled-bc} \end{eqnarray} %% where this global equation is solved for ${p^*},{\mathbf{v}^*}$, and from which ${\overline{\mathbf{v}^*}}$ can be computed. Here $\mathbf{{K}^*}$ is computed as in the corresponding {\em obstacle problem} for which the blocks are impermeable. Consider first any choice of $\mathbf{\Pi}$ of $\mathbf{\Pi}_0,\mathbf{\Pi}_1$. Consequences of Corollary~\ref{cor-ell} are that, regardless of the choice of $\mathbf{\Pi}$, the global problem is decoupled from the local problem on $\Omega_{is}$, since ${\mathbf{v}_s^*}_i$ can be written using $K_s (\mathbf{\Pi}_0 \nabla {p^*})_i$; see details below. Note that ${\mathbf{v}^*}$ is not necessarily divergence-free but ${\overline{\mathbf{v}^*}}$ is. Next we set up the upscaled version of the transport part of the system. This follows by $\nabla \cdot {\overline{\mathbf{v}^*}} =0$ and by what was said in Section~\ref{sec-newmodel-classical}. The coefficients ${\phi}^*, \mathbf{D^*}$ are defined as in the obstacle problem. % The model is as follows %% \begin{eqnarray} {\phi}^* \frac{\partial c^*}{\partial t} + {q^*}(\mathbf{x},t) - \nabla \cdot ( \mathbf{D^*} \nabla c^* - {\overline{\mathbf{v}^*}} c) &=& 0, \;\; \mathbf{x} \in \Omega, \label{eq-par-ad-global-pde} \\ % \phi_{s} \frac{\partial {c_s^*}_{i}}{\partial t} - \nabla \cdot ( \mathbf{D}_i \nabla {c_s^*}_{i} -{\mathbf{v}_s^*}_i {c_s^*}_{i}) &=& 0, \;\; \mathbf{y} \in \Omega_{is}, \label{eq-par-ad-block} \\ % {c_s^*}_i \vert_{\Gamma_i} &=& (\mathbf{\Pi} c^*)_i. \label{eq-par-ad-bc} \end{eqnarray} \end{subequations} % It remains to make a selection of $\mathbf{\Pi}$ simultaneously in the flow and in the transport parts. %%%%%%%%% \subsubsection{Case of constant approximations} Notice that with the choice $\mathbf{\Pi}=\mathbf{\Pi}_0$ in the flow equations \eqref{eq-par-flow-adv-dif}, the advection terms on $\Omega_{is}$ drop out by virtue of ${\mathbf{v}_s^*}_i(\mathbf{x}) = 0$ as in Corollary~\ref{cor-ell}. Then ${\overline{\mathbf{v}^*}}={\mathbf{v}^*}$ and is divergence-free. Next, by ${\mathbf{v}_s^*}_i(\mathbf{x}) = 0$ and \eqref{eq-dispersion}, neither advection nor dispersion effects can be captured. We have the self-adjoint discrete upscaled transport system % solved for $c^*(\mathbf{x},t),c^*_i(\mathbf{x},t)$: %% \begin{subequations} \label{eq-tra-ups}. \begin{eqnarray} % {\phi}^* \frac{\partial c^*}{\partial t} + {q^*}(\mathbf{x},t) - \nabla \cdot ( \mathbf{D^*} \nabla c^{*} - {\mathbf{v}^*} c^*) &=& 0, \;\; \mathbf{x} \in \Omega, \label{eq-tra-ups-global-pde} \\ % \phi_{i} \frac{\partial {c_s^*}_{i}}{\partial t} - \nabla \cdot ( \mathbf{D}_{i} \nabla {c_s^*}_{i} ) &=& 0, \;\; \mathbf{x} \in \Omega_{is} \label{eq-tra-ups-block} \\ % {c_s^*}_i \vert_{\Gamma_i} &=& (\mathbf{\Pi}_0 c^*)_i \label{eq-tra-ups-bc} \end{eqnarray} %% where the memory term is given by %% \begin{eqnarray} {q^*}(\mathbf{x},t) &=& \Pi'_0(\ivec{\mathbf{D}_i \nabla {c_s^*}_{i}(\mathbf{s},t) \cdot \mathbf{n}}). \label{eq-tra-ups-flux} \end{eqnarray} \end{subequations} %% In this model one captures tailing effects due to diffusion at disparate time scales but not to advection or dispersion. Advection terms are lost in the cell problem; they only appear in the global problem. Such a model was considered in \cite{Jaffre02} and in the context of thermal flow in \cite{P94,P95b}, also in \cite{Hagg04,wood.cmg}. It is interesting to note that the term ${q^*}$ could play the role of a regularizing term for large P\'eclet numbers in the global equation even though it is hard to see how this could be consistent with the absence of advection in \eqref{eq-tra-ups-block}. %%%%%%%%%%%% \subsubsection{Case of affine approximations} Now we use $\mathbf{\Pi}_1$ in both the flow part and in transport part; this will keep the advective velocities ${\mathbf{v}_s^*}_i$ from vanishing, The effect of $\mathbf{\Pi}_1$ on the elliptic part \eqref{eq-par-flow-adv-dif} was explained in Corollary~\ref{cor-ell}. We can rewrite the system \eqref{eq-par-flow-adv-dif} using a coefficient $\overline{\mathbf{{K}^*}} = \mathbf{{K}^*} + \theta_s \mathbf{K}_s$ %% \begin{subequations} \label{eq-par-flow-adv-dif-1} \begin{eqnarray} \nabla \cdot {\overline{\mathbf{v}^*}} &=& 0, \; \; \mathbf{x} \in \Omega \label{eq-par-f-global-1} \\ {\overline{\mathbf{v}^*}} &=& - \overline{\mathbf{{K}^*}} \nabla {p^*}, \label{eq-par-f-global-pde-1} %% \\ \nabla \cdot {\mathbf{v}_s^*}_i &=& 0, \; \; \mathbf{y} \in \Omega_{is}, \; i=1,\ldots \in N_{\rm incl} \\ {\mathbf{v}_s^*}_i &=& - \mathbf{K}_s \nabla {p_s^*}_i, \;\; \mathbf{y} \in \Omega_{is} \label{eq-flow2-upscaled-loc-1} \\ % {p_s^*}_i \vert_{\partial \Omega_{0s}} &=& (\mathbf{\Pi}_1({p^*}))_i %% \label{eq-flow2-upscaled-bc-1} \end{eqnarray} \end{subequations} %% Note that ${\overline{\mathbf{v}^*}}$ is divergence-free. Also, we find that ${\mathbf{v}_s^*}_i= - \mathbf{K}_s \nabla {p_s^*}_i$, is constant on each $\Omega_{is}$. In fact it is given by % \begin{eqnarray} \label{eq-loc-vel} {\mathbf{v}_s^*}_i= - \mathbf{K}_s (\mathbf{\Pi}_0 \nabla {p^*})_i = {\mathbf{K}_s}({\overline{\mathbf{{K}^*}}})^{-1} (\mathbf{\Pi}_0 {\overline{\mathbf{v}^*}})_i \end{eqnarray} % In summary, we solve \eqref{eq-par-f-global-1}, \eqref{eq-par-f-global-pde-1} and then calculate the local velocity by \eqref{eq-loc-vel}. We can now compute $\mathbf{D}_i$ from \eqref{eq-dispersion} using the (constant in $\mathbf{y}$) value of ${\mathbf{v}_s^*}_i$. In particular, $\mathbf{D}_i$ may be non-isotropic and non-diagonal. Also, it is expected that in general, ${\mathbf{v}_s^*}_i$ will vary with $i$. However, in the case considered in \cite{Hagg04}, the values ${\mathbf{v}^*}$ are essentially constant with $\mathbf{x}$ and, hence, ${\mathbf{v}_s^*}_i,\mathbf{D}_i$ do not vary with $i$. The upscaled transport system with $\mathbf{\Pi}_1$ follows as in Corollary~\ref{cor-par-pi1} %% \begin{subequations} \label{eq-tra2-ups} \begin{eqnarray} % {\phi}^* \frac{\partial c^*}{\partial t} + {q^*}(\mathbf{x},t) - \nabla \cdot ( \mathbf{D^*} \nabla c^* - {\overline{\mathbf{v}^*}} c^*) &=& 0, \;\; \mathbf{x} \in \Omega, \label{eq-tra2-ups-global-pde} \\ % \phi_{i} \frac{\partial {c_s^*}_{i}}{\partial t} - \nabla \cdot ( \mathbf{D}_{i} \nabla {c_s^*}_{i} -{\mathbf{v}_s^*}_i {c_s^*}_{i}) &=& 0, \;\; \mathbf{x} \in \Omega_{i}, \label{eq-tra2-ups-loc} \\ % {c_s^*}_i \vert_{\Gamma_i} &=& \Pi_1( c^*(\mathbf{x},t)), \label{eq-tra2-ups-bc} \end{eqnarray} %% with the memory term \begin{eqnarray} % {q^*}(\mathbf{x},t) = \mathbf{\Pi}'_1( \ivec{\mathbf{D}_i \nabla {c_s^*}_{i}(\mathbf{s},t) - {\mathbf{v}_s^*}_i(\mathbf{x}) {c_s^*}_i(\mathbf{s},t)) \cdot \mathbf{n}}, \; \mathbf{x} \in \Omega\,. \label{eq-tra2-ups-flux} \end{eqnarray} \end{subequations} %% It was observed in Corollary~\ref{cor-ell} and \cite{Arb89r} that across each block boundary $\int_{\Gamma_i} {\mathbf{v}_s^*}_i \cdot \mathbf{n} =0$. Hence, parts of the advective flux $\int_{\Gamma_i} (\mathbf{\Pi}_0 {c_s^*})_i {\mathbf{v}_s^*}_i \cdot \mathbf{n} $ vanish. However, not all advective contributions to ${q^*}$ are zero, unless, as in \cite{spagnuolo01}, the flow equations are upscaled with different operators than transport. %%%%%%%%%% \subsection{Convolution form of \eqref{eq-tra-ups} and \eqref{eq-tra2-ups}} \label{sec-realmodel} Now we rewrite \eqref{eq-tra-ups} and \eqref{eq-tra2-ups} using the representations derived in Section~\ref{sec-convolution}. This allows to partially decouple the system that is, to see the global transport equations \eqref{eq-tra-ups-global-pde}, \eqref{eq-tra2-ups-global-pde} written in terms of their global unknowns $c^*$ only. We assume for simplicity that there is an initial equilibrium in the system so that a counterpart of \eqref{eq-equilibrium} holds. We also assume that the kernels defined in \eqref{eq-kernel00}, \eqref{eq-kernelk0}, \eqref{eq-kernel0k}, \eqref{eq-kernelkk} do not vary with $i$. First we rewrite \eqref{eq-tra-ups-global-pde}. It is not difficult to see, following the proof of Proposition~\ref{prop-conv-0} and noting that ${\mathbf{v}_s^*}_i =0$, that in order to rewrite \eqref{eq-tra-ups} in a partially decoupled form, we merely need the definitions leading to the standard double porosity model \eqref{eq-par-modelf-global-pde} to which an advective term $\nabla \cdot ({\mathbf{v}^*} c^*)$ is added. Predictably, it has the form %%% \begin{multline} {\phi}^* \frac{\partial c^*}{\partial t} + \sum_i \hat{\chi}_i(\mathbf{x}) \mathcal{T}^{00}(t) \ast \ddt{(\mathbf{\Pi}_0 c^*)_i} - \nabla \cdot ( \mathbf{D^*} \nabla c^* - {\mathbf{v}^*} c^*) = 0, \;\; \mathbf{x} \in \Omega. \label{eq-tra-convolution} \end{multline} %% Next we handle \eqref{eq-tra2-ups} which is of major interest in this paper, as it preserves the advection and dispersion effects. We assume again initial equilibrium so that a counterpart of \eqref{eq-equilibrium} holds. Immediately we see that the generic solutions $r_k, k=0,1,2$ to problems \eqref{eq-def-r0}, \eqref{eq-def-rk} do not account for advection and cannot be used directly. Also, in general ${\mathbf{v}_s^*}_i$ will change with $i$. Hence we modify \eqref{eq-def-r0} appropriately as follows, for every $i=1,\ldots N_{\rm incl}$ % \begin{eqnarray} \label{eq-def-r0-adv} \left\{ \begin{array}{lll} \phi_s \frac{\partial r^0_i}{\partial t} - \nabla \cdot ( \mathbf{D}_i \nabla r^0_i - {\mathbf{v}_s^*}_i r^0_i) &= 0, \;\; &\mathbf{y} \in \Omega_{is} \,, \\ r^0_i(\mathbf{y},0) &= 0, \;\; &\mathbf{y} \in \Omega_{is} \,, \\ r^0_i(\mathbf{y},t) &= 1, \;\; &\mathbf{y} \in \Gamma_i \,. \end{array} \right. \end{eqnarray} %% We also modify \eqref{eq-def-rk} analogously, for $i=1,\ldots N_{\rm incl},k=1,2$ % \begin{eqnarray} \label{eq-def-rk-adv} \left\{ \begin{array}{lll} \phi_s \frac{\partial r^k_i}{\partial t} - \nabla \cdot ( \mathbf{D}_i \nabla r^k_i - {\mathbf{v}_s^*}_i r^0_i) &= 0, \;\; &\mathbf{y} \in \Omega_{is} \,, \\ r^k_i(\mathbf{y},0) &= 0, \;\; &\mathbf{y} \in \Omega_{is} \,, \\ r^k_i(\mathbf{y},t) &= (\mathbf{y}-\mathbf{x}^c_0)_k, \;\; &\mathbf{y} \in \Gamma_i \,. \end{array} \right. \end{eqnarray} %% Now we propose to use the definitions \eqref{eq-kernel00}, \eqref{eq-kernelk0}, \eqref{eq-kernel0k} where $r^0,r^1,r^2$ are defined by \eqref{eq-def-r0-adv}, \eqref{eq-def-rk-adv}, and allow for their variability with $i$. Also, we modify \eqref{eq-kernelkk} to include the total advective-diffusive flux $\mathbf{D}_i \nabla r^k_i - {\mathbf{v}_s^*}_i r^k_i$. This is done as follows: %% \begin{eqnarray} \label{eq-kernelkk-adv} \mathbf{S}_i^k(t) \equiv \frac{1}{\abs{\Omega_i}} \int_{\Omega_{is}} \left(\mathbf{D}_i \nabla r^k_i (\mathbf{y},t) - {\mathbf{v}_s^*}_i r^k_i(\mathbf{y},t)\right)\,dA. \end{eqnarray} %% Finally we follow calculations similar to those in Section~\ref{sec-convolution-secondary} to obtain the final result, the global upscaled discrete double-porosity model. Its structure differs from the one derived in Proposition \ref{prop-pi1} and \eqref{eq-gl-conv-pde} by the presence of the advection term $\nabla \cdot( {\mathbf{v}^*} c^*)$ and by the dependence of kernels on $i$. \begin{proposition} \label{prop-pi1-tra} The global transport equation of the discrete upscaled double-porosity model using affine approximations is given by %% \begin{multline} \label{eq-tra2-conv-pde} % {\phi}^* \dtt{c^*(\mathbf{x},t)} + \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} \mathcal{T}^{00}_i(\cdot) \ast \ddt{(\mathbf{\Pi}_0 c^*(,\cdot))_i} \\ + \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} \sum_{k=1}^2 \mathcal{T}^{k0}_i(\cdot) \ast \ddt{ (\mathbf{\Pi}_0 \partial_kc^*(\cdot))_i} \\ - \nabla \cdot \left( \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} [ \mathcal{T}^{01}_i(\cdot ),\mathcal{T}^{02}_i(\cdot)]^T \ast \ddt{ (\mathbf{\Pi}_0 c^*(\cdot))_i}\right. \\ +\left. \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} \sum_{k=1}^2 [ \mathcal{T}^{k1}_i(\cdot ),\mathcal{T}^{k2}_i(\cdot)]^T \ast \ddt{ (\mathbf{\Pi}_0 \partial_k c^*(\cdot))_i}\right. \\ \left. + \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} [S^{01}_i(\cdot),S^{02}_i(\cdot)]^T \ast \ddt{ (\mathbf{\Pi}_0 c^*(\cdot))_i } \right. \\ + \left. \sum_{i=1}^{N_{\rm incl}} {\hat{\chi}_i(\mathbf{x})} \sum_{k=0}^2 [S^{k1}_i(\cdot),S^{k2}_i(\cdot)]^T \ast \ddt{ (\mathbf{\Pi}_0\partial_k c^*(\cdot))_i } \right) \\ - \nabla \cdot \left( \mathbf{D^*} \nabla c^*(\mathbf{x},t) - {\mathbf{v}^*} c^*\right) = 0, \mathbf{x} \in \Omega,\ t > 0. \end{multline} %% \end{proposition} % This is a central result of this work. We conclude with an analogue of Remark~\ref{rem-nonpicky}. \begin{remark} \label{rem-nonpicky-tra} \rm Let $N_{\rm incl}$ be large and let us supress the dependence of the kernels on $i$. The formal limit of \eqref{eq-tra2-conv-pde} in the sense pursued in Remark~\ref{rem-nonpicky} is % \begin{multline} {\phi}^* \dtt{c^*} + \mathcal{T}^{00} \ast \dtt{c^*} + \mathbf{M} \ast \nabla \dtt{c^*} -\nabla \cdot \left( \mathcal{M} \ast \nabla \dtt{c^*} \right) - \nabla \cdot \left( \mathbf{D^*} \nabla c^* - {\mathbf{v}^*} c^* \right) = 0 \end{multline} % where $\mathbf{M}, \mathcal{M}$ have the meaning defined in Remark~\ref{rem-nonpicky}. \end{remark} The natural question that arises is one of quantitative significance of the terms associated with $\mathcal{T}^{00}, \mathbf{M}, \mathcal{M}$. Recall Remark~\ref{rem-symmetry} for the symmetric case. Numerical results and further discussion of these issues will be presented elsewhere. Further work includes construction of an $\epsilon-$model to display the upscaled model as a limit by homogenization rather than as a limit of the discrete upscaled model presented here. \begin{thebibliography}{00} \bibitem{Adams} Robert~A. Adams. \newblock {\em Sobolev spaces}. \newblock Academic Press [A subsidiary of Harcourt Brace Jovanovich, Publishers], New York-London, 1975. \newblock Pure and Applied Mathematics, Vol. 65. \bibitem{Jaffre02} Clarisse Alboin, J{\'e}r{\^o}me Jaffr{\'e}, Patrick Joly, Jean~E. Roberts, and Christophe Serres. \newblock A comparison of methods for calculating the matrix block source term in a double porosity model for contaminant transport. \newblock {\em Comput. Geosci.}, 6(3-4):523--543, 2002. \newblock Locally conservative numerical methods for flow in porous media. \bibitem{Allaire92} Gr{\'e}goire Allaire. \newblock Homogenization and two-scale convergence. \newblock {\em SIAM J. Math. Anal.}, 23(6):1482--1518, 1992. \bibitem{AGP05} B.~Amaziane, M.~Goncharenko, and L.~Pankratov. \newblock Homogenization of a degenerate triple porosity model with thin fissures. \newblock {\em European J. Appl. Math.}, 16(3):335--359, 2005. \bibitem{Arbogast} T.~Arbogast. \newblock The existence of weak solutions to single porosity and simple dual-porosity models of two-phase incompressible flow. \newblock {\em Nonlinear Analysis, Theory, Methods and Applications}, 19:1009--1031, 1992. \bibitem{Arb88} Todd Arbogast. \newblock The double porosity model for single phase flow in naturally fractured reservoirs. \newblock In {\em Numerical simulation in oil recovery (Minneapolis, Minn., 1986)}, volume~11 of {\em IMA Vol. Math. Appl.}, pages 23--45. Springer, New York, 1988. \bibitem{Arb89} Todd Arbogast. \newblock Analysis of the simulation of single phase flow through a naturally fractured reservoir. \newblock {\em SIAM J. Numer. Anal.}, 26(1):12--29, 1989. \bibitem{Arb89r} Todd Arbogast. \newblock On the simulation of incompressible, miscible displacement in a naturally fractured petroleum reservoir. \newblock {\em RAIRO Mod\'el. Math. Anal. Num\'er.}, 23(1):5--51, 1989. \bibitem{Arbogast97} Todd Arbogast. \newblock Computational aspects of dual-porosity models. \newblock In U.~Hornung, editor, {\em Homogenization and porous media}, volume~6 of {\em Interdiscip. Appl. Math.}, pages 203--215. Springer, New York, 1997. \bibitem{ADH} Todd Arbogast, Jim Douglas, Jr., and Ulrich Hornung. \newblock Derivation of the double porosity model of single phase flow via homogenization theory. \newblock {\em SIAM J. Math. Anal.}, 21(4):823--836, 1990. \bibitem{BZK60} G.~I. Barenblatt, I.~P. Zheltov, and I.~N. Kochina. \newblock Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks (strata). \newblock {\em J. Appl. Math. Mech.}, 24:1286--1303, 1960. \bibitem{Bear72} J.~Bear. \newblock {\em Dynamics of fluids in porous media}. \newblock Elsevier, New York, 1972. \bibitem{BLP} Alain Bensoussan, Jacques-Louis Lions, and George Papanicolaou. \newblock {\em Asymptotic analysis for periodic structures}, volume~5 of {\em Studies in Mathematics and its Applications}. \newblock North-Holland Publishing Co., Amsterdam, 1978. \bibitem{Bourgeat} A.~Bourgeat. \newblock Homogenized behavior of two-phase flows in naturally fractured reservoirs with uniform fractures distribution. \newblock {\em Comput. Meth. Appl. Mech. Eng.}, 47:205--216, 1984. \bibitem{Bourgeat97} Alain Bourgeat. \newblock Two-phase flow. \newblock In U.~Hornung, editor, {\em Homogenization and porous media}, volume~6 of {\em Interdiscip. Appl. Math.}, pages 97--187. Springer, New York, 1997. \bibitem{BP98} Alain Bourgeat and Mikhail Panfilov. \newblock Effective two-phase flow through highly heterogeneous porous media: capillary nonequilibrium effects. \newblock {\em Comput. Geosci.}, 2(3):191--215, 1998. \bibitem{CarslawJaeger} H.~S. Carslaw and J.~C. Jaeger. \newblock {\em Conduction of heat in solids}. \newblock Oxford Science Publications. The Clarendon Press Oxford University Press, New York, second edition, 1988. \bibitem{CiorPaulin79} Do{\"{\i}}na Cioranescu and Jeannine Saint~Jean Paulin. \newblock Homogenization in open sets with holes. \newblock {\em J. Math. Anal. Appl.}, 71(2):590--607, 1979. \bibitem{ClarkShow94} G.~W. Clark and R.~E. Showalter. \newblock Fluid flow in a layered medium. \newblock {\em Quart. Appl. Math.}, 52(4):777--795, 1994. \bibitem{CookShow95} John~D. Cook and R.~E. Showalter. \newblock Microstructure diffusion models with secondary flux. \newblock {\em J. Math. Anal. Appl.}, 189(3):731--756, 1995. \bibitem{Cushman93} J.~H. Cushman. \newblock Nonlocal dispersion in media with continuously evolving scales of heterogeneity. \newblock {\em Transp. Porous Media, 13}, 13:123--138, 1993. \bibitem{DPS97} J.~Douglas, Jr., M.~Peszy{\'n}ska, and R.~E. Showalter. \newblock Single phase flow in partially fissured media. \newblock {\em Transp. Porous Media}, 28:285--306, 1997. \bibitem{DougArb} J.~Jr. Douglas and T.~Arbogast. \newblock Dual-porosity models for flow in naturally fractured reservoirs. \newblock In J.~H. Cushman, editor, {\em Dynamics of Fluids in Hierarchical Porous Media}, pages 177--221. Academic Press, 1990. \bibitem{spagnuolo01} Jim Douglas, Jr. and Anna~M. Spagnuolo. \newblock The transport of nuclear contamination in fractured porous media. \newblock {\em J. Korean Math. Soc.}, 38(4):723--761, 2001. \newblock Mathematics in the new millennium (Seoul, 2000). \bibitem{ELW00} R.~E. Ewing, Y.~Lin, and J.~Wang. \newblock A numerical approximation of nonfickian flows with mixing length growth in porous media. \newblock {\em Acta Math. Univ. Comenian. (N.S.)}, 70(1):75--84, 2000. \bibitem{Ewing83} Richard~E. Ewing. \newblock Problems arising in the modeling of processes for hydrocarbon recovery. \newblock In Richard~E. Ewing, editor, {\em The mathematics of reservoir simulation}, volume~1 of {\em Frontiers in Applied Mathematics}, pages 3--34. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1983. \bibitem{wood.cmg} F.~Golfier, M.~Quintard, F.~Cherblanc, Harvey C.F., B.A. Zinn, and B.D. Wood. \newblock Comparison of theory and experiment for a two-region solute transport model. \newblock submitted. \bibitem{GoltzRoberts87} Mark~N. Goltz and Paul Roberts. \newblock Using the method of moments to analyze three-dimensional diffusion-limited solute transport from temporal and spatial perspectives. \newblock {\em Water Res. Research}, 23(8):1575--1585, 1987. \bibitem{roy-tailing} M.N. Gooseff, S.~M. Wondzell, R.~Haggerty, and J.~Anderson. \newblock Comparing transient storage modeling and residence time distribution (rtd) analysis in geomorphically varied reaches in the {L}ookout {C}reek basin, {O}regon, {USA}. \newblock {\em Advances in Water Resources}, 26:925--937, 2003. \bibitem{Hagg01} R.~Haggerty, S.W. Fleming, L.C. Meigs, and S.A. McKenna. \newblock Tracer tests in a fractured dolomite 2. analysis of mass transfer in single-well injection-withdrawal tests. \newblock {\em Water Resources Research}, 37(5):1129--1142, 2001. \bibitem{Hagg95} R.~Haggerty and S.M. Gorelick. \newblock Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. \newblock {\em Water Resources Research}, 31(10):2383--2400, 1995. \bibitem{Hagg00} R.~Haggerty, S.A. McKenna, and L.C. Meigs. \newblock On the late-time behavior of tracer test breakthrough curves. \newblock {\em Water Resources Research}, 36(12):3467--3479, 2000. \bibitem{roy-hyporheic} R.~Haggerty, S.~M. Wondzell, and M.~A. Johnson. \newblock Power-law residence time distribution in the hyporheic zone of a 2nd-order mountain stream. \newblock 29(13):18--1--18--4, 2002. \bibitem{Horn97} Ulrich Hornung, editor. \newblock {\em Homogenization and porous media}, volume~6 of {\em Interdisciplinary Applied Mathematics}. \newblock Springer-Verlag, New York, 1997. \bibitem{HornSho90} Ulrich Hornung and Ralph~E. Showalter. \newblock Diffusion models for fractured media. \newblock {\em J. Math. Anal. Appl.}, 147(1):69--80, 1990. \bibitem{JKO} V.~V. Jikov, S.~M. Kozlov, and O.~A. Ole{i}nik. \newblock {\em Homogenization of differential operators and integral functionals}. \newblock Springer-Verlag, Berlin, 1994. \bibitem{McleanThomee} W.~McLean and V.~Thom{\'e}e. \newblock Numerical solution of an evolution equation with a positive-type memory term. \newblock {\em J. Austral. Math. Soc. Ser. B}, 35(1):23--70, 1993. \bibitem{MclTho} W.~McLean and V.~Thom{\'e}e. \newblock Asymptotic behaviour of numerical solutions of an evolution equation with memory. \newblock {\em Asymptot. Anal.}, 14(3):257--276, 1997. \bibitem{Panfilov} Mikhail Panfilov. \newblock {\em Macroscale Models of Flow through Highly Heterogeneous Porous Media}, volume~16 of {\em Theory and Applications of Transport in Porous Media}. \newblock Kluwer Academic Publishers, Dordrecht, 2000. \bibitem{PanPiatRyb03} A.~Pankratov, A.~Piatnitskii, and V.~Rybalko. \newblock Homogenized model of reaction-diffusion in a porous medium. \newblock {\em C. R. Mecanique}, 331:253--258, 2003. \bibitem{P94} M.~Peszy\'nska. \newblock Finite element approximation of a model of nonisothermal flow through fissured media. \newblock In R.~Stenberg M.~Krizek, P.~Neittaanmaki, editor, {\em Finite Element Methods}, Lecture Notes in Pure and Applied Mathematics, pages 357--366. Marcel Dekker, 1994. \bibitem{P92} Ma{\l}gorzata Peszy\'nska. \newblock {\em Fluid flow through fissured media. Mathematical analysis and numerical approach}. \newblock PhD thesis, University of Augsburg, Augsburg, Germany, 1992. \bibitem{P95b} Ma{\l}gorzata Peszy{\'n}ska. \newblock On a model of nonisothermal flow through fissured media. \newblock {\em Differential Integral Equations}, 8(6):1497--1516, 1995. \bibitem{P95c} Ma{\l}gorzata Peszy{\'n}ska. \newblock Finite element approximation of diffusion equations with convolution terms. \newblock {\em Math. Comp.}, 65(215):1019--1037, 1996. \bibitem{P96} Ma{\l}gorzata Peszy{\'n}ska. \newblock Memory effects and microscale. \newblock In M.~Peszy\'nska K.~Malanowski, Z.~Nahorski, editor, {\em Modelling and optimization of distributed parameter systems (Warsaw, 1995)}, pages 362--369. Chapman \& Hall, New York, 1996. \bibitem{PeterBohm05} Malte~A. Peter and Michael B\"ohm. \newblock Scalings in homogenisation of reaction, diffusion and interfacial exchange in a two-phase medium. \newblock {\em Proc. Equadiff}, 11:1--8, 2005. \bibitem{Russ88} T.~Russell. \newblock Formulation of a model accounting for capillary non-equilibrium in naturally fractured subsurface formations. \newblock In R.~E. Ewing and D.~Copeland, editors, {\em Proceedings of the Fourth Wyoming Enhanced Oil Recovery Symposium}, pages 103--114, 1988. \bibitem{SP} Enrique S{\'a}nchez-Palencia. \newblock {\em Nonhomogeneous media and vibration theory}, volume 127 of {\em Lecture Notes in Physics}. \newblock Springer-Verlag, Berlin, 1980. \bibitem{ShowBook} R.~E. Showalter. \newblock {\em Hilbert space methods for partial differential equations}. \newblock Pitman, London, 1977. \newblock Monographs and Studies in Mathematics, Vol. 1. \bibitem{Show-monotone} 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{ShowVis04} R.~E. Showalter and D.~B. Visarraga. \newblock Double-diffusion models from a highly-heterogeneous medium. \newblock {\em J. Math. Anal. Appl.}, 295:191--210, 2004. \bibitem{ShowWalk91a} R.~E. Showalter and N.~J. Walkington. \newblock Diffusion of fluid in a fissured medium with microstructure. \newblock {\em SIAM J. Math. Anal.}, 22(6):1702--1722, 1991. \bibitem{ShowWalk93} R.~E. Showalter and N.~J. Walkington. \newblock Elliptic systems for a medium with microstructure. \newblock In {\em Geometric analysis and nonlinear partial differential equations (Denton, TX, 1990)}, volume 144 of {\em Lecture Notes in Pure and Appl. Math.}, pages 91--104. Dekker, New York, 1993. \bibitem{Vogt} C.~Vogt. \newblock A homogenization theorem leading to a {V}olterra integro-differential equation for permeation chromatography. \newblock Technical report, SFB 123, Heidelberg, 1982. \bibitem{WarRoo63} J.~E. Warren and P.~J. Root. \newblock The behavior of naturally fractured reservoirs. \newblock {\em Soc. Petro. Eng. Jour.}, 3:245--255, 1963. \bibitem{SanCarr04} Sanchez-Villa X. and J.~Carrera. \newblock On the striking similarity between the moments of breakthrough curves for a heterogeneous medium and a homogeneous medium with a matrix diffusion term. \newblock {\em Journal of Hydrology}, 294:164--175, 2004. \bibitem{Hagg04} Brendan Zinn, Lucy~C. Meigs, Charles~F. Harvey, Roy Haggerty, Williams~J. Peplinski, and Claudius~Freiherr von Schwerin. \newblock Experimental visualization of solute transport and mass transfer processes in two-dimensional conductivity fields with connected regions of high conductivity. \newblock {\em Environ Sci. Technol.}, 38:3916--3926, 2004. \end{thebibliography} \end{document} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%