On Limiting Behavior of Contaminant Transport Models in Coupled Surface and Groundwater Flows

Abstract

There has been a surge of work on models for coupling surface-water with groundwater flows which is at its core the Stokes-Darcy problem. The resulting (Stokes-Darcy) fluid velocity is important because the flow transports contaminants. The analysis of models including the transport of contaminants has, however, focused on a quasi-static Stokes-Darcy model. Herein we consider the fully evolutionary system including contaminant transport and analyze its quasi-static limits.

1. Introduction

The Stokes-Darcy problem describes the (slow) flow of a fluid across an interface I separating a saturated porous medium ${\text{Ω}}_p\subset {\mathbb{R}}^d$ ($d=2$ or 3) and a free flowing fluid region ${\text{Ω}}_f\subset {\mathbb{R}}^d$. Such flow is important because it transports contaminants between surface and groundwater [1,2], nutrients and oxygen between capillaries and tissue [3,4], and material in industrial filtration systems [5,6]. It also arises (at higher transport velocities) in modern fuel cells, porous combustors, advanced heat exchangers, the flow of air in the lungs and in the atmospheric boundary layer over vegetation. Adding transport involves solving one additional convection-diffusion problem with the Stokes-Darcy velocity passed from a Stokes-Darcy model and thus it seems to be a simple elaboration of the model. However, adding transport introduces new difficulties and apparently is little studied, Section 1.1.

We therefore consider the equation for the concentration $c(x,t)$ of a contaminant being transported, having a source $s(x,t)$. While each application has its own specific features, a reasonable first description of this process is the forced convection equation

$$\beta c_t+\nabla ·\left(-D\nabla c+uc\right)=s(x,t)\text{in}\text{Ω}:={\text{Ω}}_f\cup {\text{Ω}}_p\cup I$$

(1)

The free flowing fluid region’s velocity, $u_f$, and pressure, p, and the porous media’s pressure head, ϕ, and velocity, $u_p$, satisfy

$$\begin{array}{ccc}\hfill u_{f,t}-\nu \Delta u_f+\nabla p& =& f_f(x,t)\text{and}\nabla ·u_f=0\text{in}{\text{Ω}}_f\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill S_0{\varphi }_t-\nabla ·\left(K\nabla \varphi \right)& =& f_p(x,t)\text{and}{u}_p=-{\beta }^{-1}K\nabla \varphi \text{in}{\text{Ω}}_p\hfill \end{array}$$

(3)

The quasi-static limit (as $S_0\to 0$) of the predicted concentration of the full model is studied herein. The transport velocity u in the concentration Equation (1) is

$$u=\left\{\begin{array}{c}{u}_f\text{in}{\text{Ω}}_f\\ u_p\text{in}{\text{Ω}}_p\end{array}.\right.$$

(4)

For the fluid flow problem various combinations of boundary conditions on the exterior boundary $\partial \text{Ω}$ are possible and generally complicate the notation without complicating the analysis. We impose homogeneous Dirichlet boundary conditions (for clarity of exposition) and the usual initial condition

$$\begin{array}{c}{u}_f=0\text{on}\partial {\text{Ω}}_f\backslash I\text{and}\varphi =0\text{on}{\text{Ω}}_p\backslash I\\ u_f(x,0)=u_f^0\left(x\right)\text{in}{\text{Ω}}_f\text{and}\varphi (x,0)={\varphi }^0\left(x\right)\text{in}{\text{Ω}}_p\end{array}$$

For the concentration we assume that

$$\begin{array}{c}c=0\text{on}{\text{Γ}}_{in}\subset \partial \text{Ω}\text{and}-D\nabla c·\widehat{n}=0\text{on}\partial \text{Ω}\backslash \overline{{\text{Γ}}_{in}}\\ \text{and}c(x,0)=c^0\left(x\right)\text{in}\text{Ω}.\end{array}$$

There are a variety of possible interface conditions studied for I that describe different types of interfaces, e.g., [7,8,9]. Let $\widehat{n}$ be the outward pointing unit normal vector on ${\text{Ω}}_f$ and ${\left\{{\widehat{\tau }}_i\right\}}_{i=1}^{d-1}$ denote an orthonormal basis of tangent vectors on I. For slow flows across I, conservation of mass, balance of normal forces and the Beavers-Joseph-Saffman condition [10,11,12], are increasingly accepted:

$$\begin{array}{c}{u}_f·\widehat{n}-u_p·\widehat{n}=0\\ g\varphi =p-\nu \widehat{n}·\left(\nabla u_f+\nabla u_f^{\top }\right)·\widehat{n}\\ -\widehat{n}·\left(\nabla u_f+\nabla u_f^{\top }\right)·{\widehat{\tau }}_i=\frac{\alpha }{\sqrt{{\widehat{\tau }}_i·K·{\widehat{\tau }}_i}}{u}_f·{\widehat{\tau }}_i,i=1,\dots ,d-1\end{array}\left.\phantom{\begin{array}{c}{u}_f·\widehat{n}-u_p·\widehat{n}=0\\ g\varphi =p-\nu \widehat{n}·\left(\nabla u_f+\nabla u_f^{\top }\right)·\widehat{n}\\ -\widehat{n}·\left(\nabla u_f+\nabla u_f^{\top }\right)·{\widehat{\tau }}_i=\frac{\alpha }{\sqrt{{\widehat{\tau }}_i·K·{\widehat{\tau }}_i}}{u}_f·{\widehat{\tau }}_i,i=1,\dots ,d-1\end{array}}\right\}\text{on}I\text{.}$$

(5)

Interface conditions on the concentration are not needed as a single domain formulation of (1.1) imposes continuity of concentration and fluxes as natural interface conditions

$$\left[c\right]=0\text{and}[\left(-D\nabla c+uc\right)·\widehat{n}]=0,\text{on}I$$

(Jump Conditions)

The parameters in the above are as follows

$$\left\{\begin{array}{cc}{S}_0=\text{specific}\text{storage}\hfill & \nu =\text{kinematic}\text{viscosity}\hfill \\ K=\text{hydraulic}\text{conductivity}\text{tensor}\text{(SPD)}\hfill & D=\text{dispersion}\text{tensor}\hfill \\ \beta =\text{volumetric}\text{porosity}\hfill & g=\text{gravitational}\text{acceleration}\hfill \\ \alpha =\text{experimentally}\text{determined}\text{coefficient}\hfill & f_{f/p},s=\text{body}\text{forces}\text{and}\text{sources}\hfill \end{array}\right.$$

Given that $S_0$ is often very small, most of the algorithmic advances have been for the case $S_0=0$ and the concentration the primary variable of interest. The question of convergence of the concentration of the full model to that predicted by the quasi-static model as $S_0$ $\to 0$ is of significant interest and studied herein. In Theorem 1 we show that $c\to c^{QS}$ as $S_0\to 0$. This extends the analysis in [13] from the Stokes-Darcy problem to the concentration predicted by the Stokes-Darcy-Transport coupling.

The full model presents several computational and analytical difficulties (addressed herein) that are explained next. The first is an active nonlinearity in the transport problem. Taking the $L^2$ inner product of the transport equation with $c(x,t)$ and performing the standard estimates for $c(x,t)$ gives

$$\frac{1}{2}\frac{d}{dt}{\int }_{\text{Ω}}\beta c^{2}dx+{\int }_{\text{Ω}}D{\left|\nabla c\right|}^{2}dx+\frac{1}{2}{\int }_{\text{Ω}}(\nabla ·u)c^{2}dx={\int }_{\text{Ω}}scdx$$

The key term involves $\nabla ·u$ which, in the quasi-static ($S_0=0$) case, is a known function (${\beta }^{-1}{f}_p$) and, in the fully evolutionary case, is

$$\begin{array}{c}\hfill \nabla ·u=\left\{\begin{array}{c}0\text{in}{\text{Ω}}_f,\hfill \\ {\beta }^{-1}\left(-S_0\frac{\partial \varphi }{\partial t}+f_p\right)\text{in}{\text{Ω}}_p\hfill \end{array}\right.\end{array}$$

(6)

Thus, when $S_0=0$ the (nonlinear) transport term acts in the a priori estimates, stability and convergence analysis in a simpler manner than when $S_0\ne 0$.

The second issue is the multitude of small parameters in the full problem. For example, when $S_0=0$ the small parameter $K_{min}>0$ (the minimum eigenvalue of the hydraulic conductivity tensor, K) in Equation (3) can be eliminated by re-scaling $f_p.$ When $S_0\ne 0$ the small parameters in the porous media equation are active. The transport Equation (1) is also complicated by small parameters in many applications. In the simplest case where this issue occurs, it reduces to a singularly perturbed convection diffusion equation with no control on $\nabla ·u$, a problem for which methods are comparatively less well developed.

1.1. Related Work

Porous media transport and transport in a freely flowing fluid describe different physical processes with different variables, time scales, flow rates and uncertainties. There has been an intense effort at developing algorithms that use the subdomain/sub-physics codes to maximum effect to solve the coupled problem, e.g., domain decomposition methods for the equilibrium problem [14,15,16,17,18] and partitioned methods for the evolutionary problem [3,8,19,20,21]. The analytical needs to support reliability of the resulting predictions have also spurred analytical study of the coupled model. Presented in [21,22,23,24] are analyses for the coupled Stokes-Darcy-Transport problem where the velocity u in Equation (1) is modeled as either that from a fully steady Stokes-Darcy flow, or from a quasi-static coupled Stokes-Darcy flow (i.e., Equations (2) and (3) with $S_0=0$). In these the quasi-static Stokes-Darcy problem is typically solved by a domain decomposition procedure and a single domain transport problem is solved. To our knowledge, while there is for example a journal dedicated to “Transport in Porous Media”, there has been little progress on the numerical analysis of methods for full uncoupling of ( $u_p,u_f,c$) of the fully evolutionary ($S_0\ne 0$) problem.

2. Preliminaries

Let the $L^2$ norms and inner products over ${\text{Ω}}_{p/f}$ and I be denoted respectively by ${\parallel ·\parallel }_{p/f/I},{(·,·)}_{p/f/I}.$ Recall that $\text{Ω}={\text{Ω}}_p\cup {\text{Ω}}_f\cup I$; the $L^2$ norm and inner product over Ω will be denoted by $\parallel ·\parallel ,(·,·)$ (without subscripts). Let $\mathcal{D}\subset \text{Ω}$ be a regular bounded open set. With the norm ${\parallel ·\parallel }_{L^r(0,T;L^p\left(\mathcal{D}\right))}$ defined by ${\parallel f\parallel }_{L^r(0,T;L^p\left(\mathcal{D}\right))}:={\parallel \parallel f\left(t\right)\parallel }_{L^p\left(\mathcal{D}\right)}{\parallel }_{L^r(0,T)}$, we have that $L^r(0,T;L^p\left(\mathcal{D}\right)):={\{f:\parallel f\parallel }_{L^r(0,T;L^p\left(\mathcal{D}\right))}<\infty \}$.

We recall that by the Gagliardo-Nirenberg inequality [25] we have

$$\begin{array}{c}\hfill {\parallel \phi \parallel }_{L^4\left(\mathcal{D}\right)}\le C\left\{\begin{array}{cc}{\parallel \phi \parallel }_{L^2\left(\mathcal{D}\right)}^{1/2}{\parallel \phi \parallel }_{H^1\left(\mathcal{D}\right)}^{1/2}\hfill & \text{in }2d,\hfill \\ {\parallel \phi \parallel }_{L^2\left(\mathcal{D}\right)}^{1/4}{\parallel \phi \parallel }_{H^1\left(\mathcal{D}\right)}^{3/4}\hfill & \text{in }3d,\hfill \end{array}\right.\forall \phi \in H^1\left(\mathcal{D}\right)\end{array}$$

(7)

We also recall that by Remark 1.1 in [26] we have

$$\begin{array}{c}\hfill {\parallel \phi \parallel }_{L^2\left(\mathcal{D}\right)}\le C\left(\mathcal{D}\right)\left(\gamma \left(u\right)+{\parallel \nabla \phi \parallel }_{L^2\left(\mathcal{D}\right)}\right),\forall \phi \in H^1\left(\mathcal{D}\right)\end{array}$$

(8)

where $\gamma \left(u\right)$ is a seminorm, continuous on $L^2\left(\mathcal{D}\right)$, which is a norm on constants. Let Γ be a portion of $\partial \mathcal{D}$ with $meas\left(\text{Γ}\right)>0$, and assume that φ has zero trace on $\text{Γ}\subset \partial \mathcal{D}$. Then choosing $\gamma \left(u\right)={\parallel u\parallel }_{L^2\left(\text{Γ}\right)}$ we obtain from Equation (8) that the following Poincaré-Friedrichs type inequality holds on $ker\left(\text{Γ}\right)=\{\psi \in H^1\left(\mathcal{D}\right);\psi {|}_{\text{Γ}}=0\}$:

$$\begin{array}{c}\hfill {\parallel \phi \parallel }_{L^2\left(\mathcal{D}\right)}\le C\left(\mathcal{D}\right){\parallel \nabla \phi \parallel }_{L^2\left(\mathcal{D}\right)},\forall \phi \in H^1\left(\mathcal{D}\right),\phi {|}_{\text{Γ}}=0\end{array}$$

(9)

From Equations (7) and (9) we derive

$$\begin{array}{c}\hfill {\parallel \phi \parallel }_{L^4\left(\mathcal{D}\right)}\le C\left\{\begin{array}{cc}{\parallel \phi \parallel }_{L^2\left(\mathcal{D}\right)}^{1/2}{\parallel \nabla \phi \parallel }_{L^2\left(\mathcal{D}\right)}^{1/2}\hfill & \text{in }2d,\hfill \\ {\parallel \phi \parallel }_{L^2\left(\mathcal{D}\right)}^{1/4}{\parallel \nabla \phi \parallel }_{L^2\left(\mathcal{D}\right)}^{3/4}\hfill & \text{in }3d,\hfill \end{array}\right.\forall \phi \in H^1\left(\mathcal{D}\right),\phi {|}_{\text{Γ}}=0\end{array}$$

(10)

If $\phi =0$ on $\partial \mathcal{D}$, then Equation (10) are just the inequalities proved by Ladyzhenskaya [27]. Denote the (assumed positive) minimum of D by

$$D_{min}=\underset{x\in \text{Ω}}{inf}D\left(x\right)>0$$

Regularity of the concentration depends on regularity of the Stokes-Darcy variables. In [13] Moraiti proved that for $0<T<\infty$ and data satisfying

$$\begin{array}{cc}\hfill f_{f,t}\in L^2(0,T;H^{-1}\left({\text{Ω}}_f\right))& ,f_{p,t}\in L^2(0,T;H^{-1}\left({\text{Ω}}_p\right))\hfill \\ \hfill u_{f,t}\left(0\right)\in L^2\left({\text{Ω}}_f\right)& ,{\varphi }_t\left(0\right)\in L^2\left({\text{Ω}}_p\right)\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill u_{f,t}\left(0\right)& :=u_{f,t}(x,0)=\underset{t\to 0^{+}}{lim}{u}_{f,t}(x,t)=\underset{t\to 0^{+}}{lim}\left(f_f(x,t)+\nu \Delta u_f(x,t)-\nabla p(x,t)\right)\hfill \\ \hfill {\varphi }_t\left(0\right)& :={\varphi }_t(x,0)=\underset{t\to 0^{+}}{lim}{\varphi }_t(x,t)=S_0^{-1}\underset{t\to 0^{+}}{lim}\left(f_p(x,t)+\nabla ·\left(K\nabla \varphi (x,t)\right)\right)\hfill \end{array}$$

the following hold uniformly in $S_0$ and will be assumed herein:

$$u_{f,t}\in L^{\infty }(0,T;L^2\left({\text{Ω}}_f\right)),\sqrt{S_0}{\varphi }_t\in L^{\infty }(0,T;L^2\left({\text{Ω}}_p\right))\text{and}\nabla {\varphi }_t\in L^2(0,T;L^2\left({\text{Ω}}_p\right))$$

(11)

Additionally we assume

$$\begin{array}{cc}& c^0\in L^2\left(\text{Ω}\right),\nabla c^0\in L^2\left(\text{Ω}\right),s\in L^2(0,T;L^2\left(\text{Ω}\right)),\nabla {\varphi }_t\left(0\right)\in L^2\left({\text{Ω}}_p\right),\text{and}\hfill \end{array}$$

(12)

$$\begin{array}{cc}& f_f\in L^{\infty }(0,T;L^2\left({\text{Ω}}_f\right)),f_p\in L^{\infty }(0,T;H^1\left({\text{Ω}}_p\right)),f_{p,t}\in L^2(0,T;L^2\left({\text{Ω}}_p\right))\hfill \end{array}$$

(13)

We note that from the Stokes Equation (2) we have

$$\begin{array}{c}\hfill \parallel \Delta u_f{\parallel }_{L^{\infty }(0,T;L^2\left({\text{Ω}}_f\right))}\le C\left(\parallel u_{f,t}{\parallel }_{L^{\infty }(0,T;L^2\left({\text{Ω}}_f\right))}^2+{\parallel f_f\parallel }_{L^{\infty }(0,T;L^2\left({\text{Ω}}_f\right))}^2\right)\end{array}$$

hence under the regularity assumptions in Equations (11) and (13) we obtain that $u_f\in L^{\infty }(0,T;H^2\left({\text{Ω}}_f\right))$. Also, using energy estimate arguments similar to [13] for the Darcy equation, we have

$$\begin{array}{c}\hfill S_0\parallel K^{1/2}\nabla {\varphi }_t{\left(t\right)\parallel }_p^2+{\int }_0^t\parallel \nabla ·\left(K\nabla {\varphi }_t\right){\parallel }_p^{2}dr\le \parallel f_{t,p}{\parallel }_{L^2(0,T;L^2\left({\text{Ω}}_p\right))}^2+S_0{\parallel K^{1/2}\nabla {\varphi }_t\left(0\right)\parallel }_p^2\end{array}$$

which under hypotheses Equations (12) and (13) gives $\sqrt{S_0}\nabla {\varphi }_t\in L^{\infty }(0,T;L^2\left({\text{Ω}}_p\right))$. Moreover, if K satisfies additional assumptions (e.g., piecewise constant or sufficiently smooth), this also implies $K\nabla \varphi \in L^{\infty }(0,T;H^2\left({\text{Ω}}_p\right))$. To summarize, in the remainder we assume that uniformly in $S_0$

$$\begin{array}{c}\hfill \begin{array}{c}{u}_f\in W^{1,\infty }([0,T];L^2\left({\text{Ω}}_f\right))\cap L^{\infty }(0,T;H^2\left({\text{Ω}}_f\right)),\hfill \\ \sqrt{S_0}\varphi \in W^{1,\infty }([0,T];H^1\left({\text{Ω}}_p\right)),K\nabla \varphi \in L^{\infty }(0,T;H^2\left({\text{Ω}}_p\right))\hfill \end{array}\end{array}$$

(14)

and we shall prove in Propositions 1 and 2 that

$$c\in \{g:g\in L^{\infty }(0,T;H^1\left(\text{Ω}\right))\cap W^{1,2}([0,T];L^2\left(\text{Ω}\right)),g{|}_{{\text{Γ}}_{in}}=0\}$$

and give estimates of $\parallel c-c^{QS}\parallel$ as $S_0\to 0$. Throughout we use C to denote a generic positive constant, whose actual value may vary from line to line in the analysis. We begin with the following a priori estimate.

Proposition 1. (The first estimate) Suppose $0<T<\infty$, the problem data for Equations (2) and (3) is such that Equations (11)–(13) hold. Then

$$c\in L^{\infty }(0,T;L^2\left(\text{Ω}\right))\text{and}\nabla c\in L^2(0,T;L^2\left(\text{Ω}\right)).$$

(15)

Proof. For the transport Equation (1), multiply by $c(x,t)$ and integrate over Ω. This gives

$$\frac{1}{2}\beta \frac{d}{dt}{\parallel c\parallel }^2+{\parallel \sqrt{D}\nabla c\parallel }^2+\frac{1}{2}{\int }_{\text{Ω}}(\nabla ·u)c^{2}dx=(s,c)$$

Since $\nabla ·u=0$ in the fluid region and $\nabla ·u=-{\beta }^{-1}\left(S_0{\varphi }_t-f_p\right)$ in the porous media region the third term is

$$\frac{1}{2}{\int }_{\text{Ω}}(\nabla ·u)c^{2}dx=-\frac{1}{2\beta }{\int }_{{\text{Ω}}_p}\left(S_0{\varphi }_t-f_p\right)c^{2}dx$$

Thus we have

$$\begin{array}{ccc}\hfill \frac{1}{2}\beta \frac{d}{dt}{\parallel c\parallel }^2+{\parallel \sqrt{D}\nabla c\parallel }^2& =& (s,c)-\frac{1}{2\beta }{\int }_{{\text{Ω}}_p}{f}_p{c}^{2}dx+\frac{S_0}{2\beta }{\int }_{{\text{Ω}}_p}{\varphi }_t{c}^{2}dx\hfill \\ & \le & \frac{1}{2}{\parallel s\parallel }^2+\frac{1}{2}{\parallel c\parallel }^2-\frac{1}{2\beta }{\int }_{{\text{Ω}}_p}{f}_p{c}^{2}dx+\frac{S_0}{2\beta }{\int }_{{\text{Ω}}_p}{\varphi }_t{c}^{2}dx\hfill \end{array}$$

The critical term is $\int {\varphi }_t{c}^{2}dx$ and estimates for this term depend on estimates for ${\varphi }_t$. Thus, by Hölder’s inequality

$$\left|\frac{S_0}{2\beta }{\int }_{{\text{Ω}}_p}{\varphi }_t{c}^{2}dx\right|\le \frac{S_0}{2\beta }\parallel {\varphi }_t{\parallel }_p{\parallel c\parallel }_{L^4\left({\text{Ω}}_p\right)}^2$$

Inequalities Equation (10) for ${\parallel c\parallel }_{L^4}^2$ imply

$$\left|\frac{S_0}{2\beta }{\int }_{{\text{Ω}}_p}{\varphi }_t{c}^{2}dx\right|\le \frac{C}{\beta }{S}_0{\parallel {\varphi }_t\parallel }_p\left\{\begin{array}{cc}{\parallel c\parallel }_p{\parallel \nabla c\parallel }_p\hfill & \text{in}2d\text{,}\hfill \\ {\parallel c\parallel }_p^{1/2}{\parallel \nabla c\parallel }_p^{3/2}\hfill & \text{in}3d\hfill \end{array}\right.$$

(16)

An analogous bound to Equation (16) holds for ${\int }_{{\text{Ω}}_p}{f}_p{c}^{2}dx$.

We consider the $2d$ and $3d$ cases separately.

The 2d case. Since ${\parallel \nabla c\parallel }_p\le D_{min}^{-1/2}{\parallel \sqrt{D}\nabla c\parallel }_p$, in $2d$ we have

$$\begin{array}{cc}& \frac{1}{2}\beta \frac{d}{dt}{\parallel c\parallel }^2+{\parallel \sqrt{D}\nabla c\parallel }^2\hfill \\ & \le \frac{1}{2}{\parallel s\parallel }^2+\frac{1}{2}{\parallel c\parallel }^2+\frac{C}{\beta }\left(\parallel f_p{\parallel }_p+S_0{\parallel {\varphi }_t\parallel }_p\right){\parallel c\parallel }_p{\parallel \nabla c\parallel }_p\hfill \\ & \le \frac{1}{2}{\parallel s\parallel }^2+\frac{1}{2}{\parallel c\parallel }^2+\frac{1}{2}\parallel \sqrt{D}{\nabla c\parallel }_p^2+D_{min}^{-1}\frac{C}{{\beta }^2}\left(\parallel f_p{\parallel }_p^2+S_0^2{\parallel {\varphi }_t\parallel }_p^2\right){\parallel c\parallel }_p^2\hfill \end{array}$$

Thus we have

$$\begin{array}{cc}& \frac{d}{dt}{\parallel c\parallel }^2+\frac{1}{\beta }\parallel \sqrt{D}{\nabla c\parallel }^2\le \frac{1}{\beta }{\parallel s\parallel }^2+\frac{1}{\beta }\left(1+D_{min}^{-1}\frac{C}{{\beta }^2}\left(\parallel f_p{\parallel }_p^2+S_0^2{\parallel {\varphi }_t\parallel }_p^2\right)\right){\parallel c\parallel }^2\hfill \end{array}$$

(17)

Proceeding as in the proof of Grönwall’s inequality, with $\mu \left(t\right)={\int }_0^t\frac{1}{\beta }\left[1+D_{min}^{-1}\frac{C}{{\beta }^2}(\parallel f_p{\left(\xi \right)\parallel }_p^2+S_0^2\parallel {\varphi }_t\left(\xi \right){\parallel }_p^2)\right]d\xi$, multiplying Equation (17) by $exp(-\mu (t\left)\right)$ and rearranging we have

$$\frac{d}{dt}\left({exp(-\mu \left(t\right))\parallel c\left(t\right)\parallel }^2\right)+\frac{1}{\beta }exp(-\mu \left(t\right))\parallel \sqrt{D}{\nabla c\parallel }^2\le \frac{1}{\beta }exp(-\mu \left(t\right)){\parallel s\parallel }^2$$

(18)

Integrating Equation (18) from 0 to t, and then multiplying through by $exp\left(\mu \right(t\left)\right)$ yields

$$\begin{array}{cc}& {\parallel c\left(t\right)\parallel }^2+\frac{1}{\beta }{\int }_0^{t}exp(\mu \left(t\right)-\mu \left(\xi \right)){\parallel \sqrt{D}\nabla c\left(\xi \right)\parallel }^{2}d\xi \hfill \end{array}$$

(19)

$$\begin{array}{cc}& \le exp\left(\mu \left(t\right)\right)\parallel c^0{\parallel }^2+\frac{1}{\beta }{\int }_0^{t}exp(\mu \left(t\right)-\mu \left(\xi \right)){\parallel s\left(\xi \right)\parallel }^{2}d\xi \hfill \end{array}$$

With the assumed regularity Equations (11) and (12), and the boundedness of $exp\left(\mu \right(T\left)\right)$, Equation (15) now follows.

The 3d case. In $3d$ we have

$$\begin{array}{cc}& \frac{1}{2}\beta \frac{d}{dt}{\parallel c\parallel }^2+{\parallel \sqrt{D}\nabla c\parallel }^2\hfill \\ & \le \frac{1}{2}{\parallel s\parallel }^2+\frac{1}{2}{\parallel c\parallel }^2+\frac{C}{\beta }\left(\parallel f_p{\parallel }_p+S_0{\parallel {\varphi }_t\parallel }_p\right){\parallel c\parallel }_p^{1/2}{\parallel \nabla c\parallel }_p^{3/2}\hfill \\ & \le \frac{1}{2}{\parallel s\parallel }^2+\frac{1}{2}{\parallel c\parallel }^2+{\parallel \sqrt{D}\nabla c\parallel }_p^{3/2}\left(D_{min}^{-3/4}\frac{C}{\beta }(\parallel f_p{\parallel }_p+S_0\parallel {\varphi }_t{\parallel }_p{)\parallel c\parallel }_p^{1/2}\right)\hfill \end{array}$$

For the last term we use $ab\le \frac{3}{4}{a}^{4/3}+\frac{1}{4}{b}^4$. This gives, after rearranging,

$$\begin{array}{c}\frac{d}{dt}{\parallel c\parallel }^2+\frac{1}{2\beta }\parallel \sqrt{D}{\nabla c\parallel }^2\le \frac{1}{\beta }{\parallel s\parallel }^2+\frac{1}{\beta }\left(1+D_{min}^{-3}\frac{C}{{\beta }^4}\left(\parallel f_p{\parallel }_p^4+S_0^4{\parallel {\varphi }_t\parallel }_p^4\right)\right){\parallel c\parallel }^2\end{array}$$

Now, proceeding as in the 2d case we obtain Equation (15).  ☐

3. Validity of the Quasi-Static Model

Let $c^{QS}(x,t)$ be the solution of Equation (1) with $S_0=0$, i.e., $u=u^{QS}$, the solution of the quasi-static Stokes-Darcy problem, where

$$\nabla ·u^{QS}=\left\{\begin{array}{c}0,\text{in}{\text{Ω}}_f,\hfill \\ \frac{1}{\beta }{f}_p,\text{in}{\text{Ω}}_p\hfill \end{array}\right.$$

(20)

Define

$$e^c(x,t):=c(x,t)-c^{QS}(x,t)\text{and}{e}^u(x,t):=u(x,t)-u^{QS}(x,t)$$

and note that $e^c(x,0)=0$, and $e^u(x,0)=0$. In Theorem 1 we show that $c\to c^{QS}$ as $S_0\to 0$. To prove convergence in $3d$, we first obtain a second a priori bound for the concentration c, given next.

Proposition 2. (The second estimate) Assuming Equations (11) and (12), we have that uniformly in $S_0$

$$\nabla c\in L^{\infty }(0,T;L^2\left(\text{Ω}\right))\text{and}{c}_t\in L^2(0,T;L^2\left(\text{Ω}\right))$$

(21)

Proof. Take the inner product of Equation (1) with $c_t$, integrate over Ω, and apply the divergence theorem. This yields:

$$\begin{array}{c}\beta \left(c_t,c_t\right)-\left(\nabla ·\left(D\nabla c\right),c_t\right)+\left(\nabla ·\left(uc\right),c_t\right)=\left(s,c_t\right)\text{and}\text{thus}\\ \beta \parallel c_t{\parallel }^2+\frac{1}{2}\frac{d}{dt}{\parallel \sqrt{D}\nabla c\parallel }^2=\left(s,c_t\right)-\left(\nabla ·\left(uc\right),c_t\right)\end{array}$$

Using Cauchy-Schwarz and Young inequalities and absorbing terms on the left-hand side, we have

$$\begin{array}{cc}\hfill \frac{\beta }{2}\parallel c_t{\parallel }^2+\frac{1}{2}\frac{d}{dt}{\parallel \sqrt{D}\nabla c\parallel }^2& \le {\beta }^{-1}\left({\parallel s\parallel }^2+{\parallel \nabla ·\left(uc\right)\parallel }^2\right)\hfill \\ & ={\beta }^{-1}\left({\parallel s\parallel }^2+{\parallel c\nabla ·u+u·\nabla c\parallel }^2\right)\hfill \\ & \le {\beta }^{-1}\left({\parallel s\parallel }^2{+2\parallel c\nabla ·u\parallel }^2+2{\parallel u·\nabla c\parallel }^2\right)\hfill \end{array}$$

We treat only the $3d$ case, because the $2d$ case follows in a similar way. Integrating over $(0,t)$, $0<t\le T$, using the Young inequality and inequality Equation (10) we get

$$\begin{array}{cc}& \parallel \sqrt{D}{\nabla c\left(t\right)\parallel }^2+\beta {\int }_0^t{\parallel c_r\left(r\right)\parallel }^{2}dr\hfill \\ & \le \parallel \sqrt{D}{\nabla c\left(0\right)\parallel }^2+2{\beta }^{-1}{\int }_0^t\left({\parallel s\left(r\right)\parallel }^2{+2\parallel c\left(r\right)\nabla ·u\left(r\right)\parallel }^2+2{\parallel u\left(r\right)·\nabla c\left(r\right)\parallel }^2\right)dr\hfill \\ & \le \parallel \sqrt{D}{\nabla c\left(0\right)\parallel }^2+2{\beta }^{-1}{\int }_0^t\left({\parallel s\left(r\right)\parallel }^2+{2\parallel c\left(r\right)\parallel }_{L^4\left(\text{Ω}\right)}^2{\parallel \nabla ·u\left(r\right)\parallel }_{L^4\left(\text{Ω}\right)}^2+{2\parallel u\left(r\right)\parallel }_{L^{\infty }\left(\text{Ω}\right)}^2{\parallel \nabla c\left(r\right)\parallel }^2\right)dr\hfill \\ & \le \parallel \sqrt{D}{\nabla c\left(0\right)\parallel }^2+2{\beta }^{-1}{\int }_0^t{\parallel s\left(r\right)\parallel }^{2}dr+4{\beta }^{-1}{\int }_0^t\left({\parallel \nabla c\left(r\right)\parallel }^2+{C\parallel c\left(r\right)\parallel }^2{\parallel \nabla ·u\left(r\right)\parallel }_{L^4\left(\text{Ω}\right)}^8\right)dr\hfill \\ & +4{\beta }^{-1}{\int }_0^t{\parallel u\left(r\right)\parallel }_{L^{\infty }\left(\text{Ω}\right)}^2{\parallel \nabla c\left(r\right)\parallel }^{2}dr\hfill \\ & \le \parallel \sqrt{D}{\nabla c\left(0\right)\parallel }^2+2{\beta }^{-1}{\parallel s\parallel }_{L^2(0,T;L^2\left(\text{Ω}\right))}^2+4{\beta }^{-1}{\parallel \nabla c\parallel }_{L^2(0,T;L^2\left(\text{Ω}\right))}^2\hfill \\ \text{(22)}\hfill & & +C{\beta }^{-1}{\int }_0^t{\parallel c\left(r\right)\parallel }^2{\parallel \nabla ·u\left(r\right)\parallel }_{L^4\left(\text{Ω}\right)}^{8}dr+4{\beta }^{-1}{\parallel u\parallel }_{L^{\infty }(0,T;L^{\infty }\left(\text{Ω}\right))}^2{\parallel \nabla c\parallel }_{L^2(0,T;L^2\left(\text{Ω}\right))}^2\hfill \end{array}$$

The second to last term in Equation (22) is treated as follows. From Equation (6) and again using Equation (10) we have

$$\begin{array}{cc}& {\int }_0^t{\parallel c\left(r\right)\parallel }^2{\parallel \nabla ·u\left(r\right)\parallel }_{L^4\left(\text{Ω}\right)}^{8}dr\le {\parallel c\parallel }_{L^{\infty }(0,T;L^2\left(\text{Ω}\right))}^2{\int }_0^t{\parallel {\beta }^{-1}(-S_0{\varphi }_t\left(r\right)+f_p\left(r\right))\parallel }_{L^4\left({\text{Ω}}_p\right)}^{8}dr\hfill \\ & \le {\beta }^{-8}{2}^7{\parallel c\parallel }_{L^{\infty }(0,T;L^2\left(\text{Ω}\right))}^2{\int }_0^t\left(S_0^8\parallel {\varphi }_t{\left(r\right)\parallel }_{L^4\left({\text{Ω}}_p\right)}^8+{\parallel f_p\left(r\right)\parallel }_{L^4\left({\text{Ω}}_p\right)}^8\right)dr\hfill \\ & \le {C\parallel c\parallel }_{L^{\infty }(0,T;L^2\left(\text{Ω}\right))}^2{\int }_0^t\left(S_0^8\parallel {\varphi }_t{\left(r\right)\parallel }_p^2\parallel \nabla {\varphi }_t{\left(r\right)\parallel }_p^6+\parallel f_p{\parallel }_p^2{\parallel \nabla f_p\parallel }_p^6\right)dr\hfill \\ & \le {C\parallel c\parallel }_{L^{\infty }(0,T;L^2\left(\text{Ω}\right))}^2\left(S_0^8\parallel {\varphi }_t{\parallel }_{L^{\infty }(0,T;L^2\left({\text{Ω}}_p\right))}^2{\int }_0^t\parallel \nabla {\varphi }_t{\left(r\right)\parallel }_p^{6}dr+\parallel f_p{\parallel }_{L^{\infty }(0,T;L^2\left({\text{Ω}}_p\right))}^2{\int }_0^t{\parallel \nabla f_p\left(r\right)\parallel }_p^{6}dr\right).\hfill \end{array}$$

For estimating the norm ${\parallel u\parallel }_{L^{\infty }(0,T;L^2\left(\text{Ω}\right))}^2$ in the last term in Equation (22) we use Equations (4) and (3), Sobolev embeddings and Equation (14)

$$\begin{array}{cc}& {\parallel u\parallel }_{L^{\infty }(0,T;L^{\infty }\left(\text{Ω}\right))}^2=\parallel u_f{\parallel }_{L^{\infty }(0,T;L^{\infty }\left({\text{Ω}}_f\right))}^2+{\parallel u_p\parallel }_{L^{\infty }(0,T;L^{\infty }\left({\text{Ω}}_p\right))}^2\hfill \\ & =\parallel u_f{\parallel }_{L^{\infty }(0,T;L^{\infty }\left({\text{Ω}}_f\right))}^2+{\beta }^{-2}{\parallel K\nabla \varphi \parallel }_{L^{\infty }(0,T;L^{\infty }\left({\text{Ω}}_p\right))}^2\hfill \\ & \le C\left(\parallel u_f{\parallel }_{L^{\infty }(0,T;H^2\left({\text{Ω}}_f\right))}^2+{\beta }^{-2}{\parallel K\nabla \varphi \parallel }_{L^{\infty }(0,T;H^2\left({\text{Ω}}_p\right))}^2\right).\hfill \end{array}$$

Finally, using Equation (12)–(15), and taking the supremum over $[0,T]$, we obtain Equation (21). ☐

We can now prove convergence of the concentration to the quasi-static approximation.

Theorem 1. (Quasi-static limit) Assume Equations (11) and (12) hold. Then for $T<\infty$

$$\begin{array}{cc}\hfill \parallel e^c{\parallel }_{L^{\infty }(0,T;L^2\left(\text{Ω}\right))}& \le C(T,data)S_0\hfill \\ \hfill \parallel \sqrt{D}\nabla e^c{\parallel }_{L^2(0,T;L^2\left(\text{Ω}\right))}& \le C(T,data)S_0\hfill \end{array}$$

Proof. Subtract the concentration equation and its quasi-static form. Next add and subtract $u^{QS}c$ in the transport term ($\nabla ·\left(uc-u^{QS}{c}^{QS}\right)$):

$$\begin{array}{cc}& \beta e_t^c-\nabla ·\left(D\nabla e^c\right)+\nabla ·\left(uc-u^{QS}{c}^{QS}\right)=0\text{and}\text{thus}\hfill \\ & \beta e_t^c-\nabla ·\left(D\nabla e^c\right)+\nabla ·\left(e^{u}c\right)+\nabla ·\left(u^{QS}{e}^c\right)=0\hfill \end{array}$$

Take the inner product with $e^c$, integrate over Ω, and apply integration by parts to obtain

$$\frac{\beta }{2}\frac{d}{dt}\parallel e^c{\parallel }^2+{\parallel \sqrt{D}\nabla e^c\parallel }^2-\left(c,e^u·\nabla e^c\right)+\left(\nabla ·\left(u^{QS}{e}^c\right),e^c\right)=0$$

Expanding and using integration by parts, we write

$$\begin{array}{cc}\hfill \left(\nabla ·\left(u^{QS}{e}^c\right),e^c\right)& =\frac{1}{2}\left(\nabla ·\left(u^{QS}{e}^c\right),e^c\right)+\frac{1}{2}\left(\nabla ·\left(u^{QS}{e}^c\right),e^c\right)\hfill \\ & =\frac{1}{2}(\nabla ·u^{QS},{\left(e^c\right)}^2)+\cancel{\frac{1}{2}\left(u^{QS},e^c\nabla e^c\right)}+\frac{1}{2}{\langle u^{QS}·\widehat{n},0\overline{){\left(e^c\right)}^2}\rangle }_{\partial \text{Ω}}^0-\cancel{\frac{1}{2}\left(u^{QS},e^c\nabla e^c\right)}\hfill \\ & =\frac{1}{2\beta }\left(f_p,{\left(e^c\right)}^2\right)\hfill \end{array}$$

Hence,

$$\frac{\beta }{2}\frac{d}{dt}\parallel e^c{\parallel }^2+{\parallel \sqrt{D}\nabla e^c\parallel }^2=\left(c,e^u·\nabla e^c\right)-\frac{1}{2\beta }\left(f_p,{\left(e^c\right)}^2\right)$$

(23)

Applying Equation (10), Poincaré-Friedrichs and Young’s inequalities we bound

$$\begin{array}{cc}\hfill \left(c,e^u·\nabla e^c\right)& \le C\left\{\begin{array}{c}\parallel \nabla e^c\parallel \parallel e^u{\parallel }^{1/2}\parallel \nabla e^u{\parallel }^{1/2}{\parallel c\parallel }^{1/2}{\parallel \nabla c\parallel }^{1/2},\text{in}\text{2d}\hfill \\ \parallel \nabla e^c\parallel \parallel e^u{\parallel }^{1/4}\parallel \nabla e^u{\parallel }^{3/4}{\parallel c\parallel }^{1/4}{\parallel \nabla c\parallel }^{3/4},\text{in}\text{3d}\hfill \end{array}\right.\hfill \\ & \le C\left\{\begin{array}{c}{D}_{min}^{-1/2}\parallel \sqrt{D}\nabla e^c\parallel \parallel \nabla e^u\parallel \parallel \nabla c\parallel ,\text{in}\text{2d}\hfill \\ D_{min}^{-1/2}\parallel \sqrt{D}\nabla e^c\parallel \parallel \nabla e^u\parallel \parallel \nabla c\parallel ,\text{in}\text{3d}\hfill \end{array}\right.\hfill \\ \text{(24)}\hfill & & \le \frac{1}{4}\parallel \sqrt{D}\nabla e^c{\parallel }^2+C\parallel \nabla e^u{\parallel }^2{\parallel \nabla c\parallel }^2\hfill \end{array}$$

Next,

$$\begin{array}{cc}\hfill \frac{1}{2\beta }\left(f_p,{\left(e^c\right)}^2\right)& \le \frac{C}{2\beta }{\parallel f_p\parallel }_p\left\{\begin{array}{c}\parallel e^c\parallel \parallel \nabla e^c\parallel ,\text{in}\text{2d}\hfill \\ \parallel e^c{\parallel }^{1/2}{\parallel \nabla e^c\parallel }^{3/2},\text{in}\text{3d}\hfill \end{array}\right.\hfill \\ \text{(25)}\hfill & & \le \left\{\begin{array}{c}\frac{1}{4}\parallel \sqrt{D}\nabla e^c{\parallel }^2+\frac{C}{{\beta }^2}{D}_{min}^{-1}\parallel f_p{\parallel }_p^2{\parallel e^c\parallel }^2,\text{in}\text{2d}\hfill \\ \frac{1}{4}\parallel \sqrt{D}\nabla e^c{\parallel }^2+\frac{C}{4{\beta }^4}{D}_{min}^{-3}\parallel f_p{\parallel }_p^4{\parallel e^c\parallel }^2,\text{in}\text{3d}\hfill \end{array}\right.\hfill \end{array}$$

We focus on the 3d case. The 2d case follows similarly. Combining Equations (23)–(25), and rearranging we have

$$\frac{d}{dt}\parallel e^c{\parallel }^2+\frac{1}{\beta }\parallel \sqrt{D}\nabla e^c{\parallel }^2\le \frac{C}{{\beta }^5}{D}_{min}^{-3}\parallel f_p{\parallel }_p^4\parallel e^c{\parallel }^2+C\parallel \nabla e^u{\parallel }^2{\parallel \nabla c\parallel }^2$$

i.e., for $\mu \left(t\right)=C{D}_{min}^{-3}/\left({\beta }^5\right){\int }_0^t{\parallel f_p\left(\xi \right)\parallel }_p^{4}d\xi$, and $\parallel e^c{\left(0\right)\parallel }^2=0$,

$$\parallel e^c{\left(t\right)\parallel }^2+\frac{1}{\beta }{\int }_0^{t}exp(\mu \left(t\right)-\mu \left(\xi \right))\parallel \sqrt{D}\nabla e^c{\left(\xi \right)\parallel }^{2}d\xi \le C{\int }_0^{t}exp(\mu \left(t\right)-\mu \left(\xi \right))\parallel \nabla e^u{\left(\xi \right)\parallel }^2{\parallel \nabla c\left(\xi \right)\parallel }^{2}d\xi$$

(26)

In [13] it is proven that, under the stated assumptions, the following hold for $e_f^u:=u_f-u_f^{QS}$, $e_p^u:=u_p-u_p^{QS}$:

$$\begin{array}{cc}\hfill \parallel e_f^u{\parallel }_{L^{\infty }(0,T;L^2\left({\text{Ω}}_f\right))}& =\mathcal{O}\left(S_0\right)\text{,}{\parallel e_p^u\parallel }_{L^{\infty }(0,T;L^2\left({\text{Ω}}_p\right))}=\mathcal{O}\left(\sqrt{S_0}\right)\hfill \\ \hfill \parallel \nabla e_f^u{\parallel }_{L^2(0,T;L^2\left({\text{Ω}}_f\right))}& =\mathcal{O}\left(S_0\right)\text{,}{\parallel \nabla e_p^u\parallel }_{L^2(0,T;L^2\left({\text{Ω}}_p\right))}=\mathcal{O}\left(S_0\right)\hfill \end{array}$$

and the analysis revealed that the convergence is sensitive in $K_{min}$, in that the constants in the convergence analysis are proportional to $1/\sqrt{K_{min}}$. Thus,

$$\parallel \nabla e^u{\parallel }_{L^2(0,T;L^2\left(\text{Ω}\right))}=\mathcal{O}\left(S_0\right)$$

(27)

With the a priori bound in Equation (15), our assumptions Equation (12), the boundedness of $exp\left(\mu \right(T\left)\right)$, taking the supremum over $[0,T]$ in Equation (26), in view of Equation (27), we obtain the first-order convergence of c to $c^{QS}$ as $S_0\to 0$, completing the proof. ☐

4. Conclusions

We conclude that the quasi-static transport model for the concentration of contaminants is justified when the specific storage parameter, $S_0$, is small when compared to the minimum eigenvalues $K_{min}$ and $D_{min}$ of the hydraulic conductivity tensor, K, and dispersion tensor, D, respectively.