2. The Stokes–Darcy equations

In this section, we show that the Stokes–Darcy system together with appropriate interface boundary conditions can be recovered from the Helmholtz minimal dissipation principle (Batchelor Reference Batchelor2000) provided we identify appropriate energy dissipation rates. We then present a non-dimensionalized version of the coupled system using the typical velocity and length of the conduit (free zone) and introduce our key small parameter—the Darcy number—together with the other physical non-dimensional parameters.

The Stokes system, the Darcy equation and the BJ interface boundary condition are used to motivate the dissipation functions but not in any other manner. Recall that different systems may enjoy the same energy dissipation rate. Indeed, Fields Medalist Terrence Tao constructed a system closely related to the three-dimensional Navier–Stokes system with the same energy dissipation rate. However, the solution to the modified system blows up in finite time while the same question for the three-dimensional Navier–Stokes equation is one of the seven open questions for a million-dollar Millennium Prize posed by the Clay Mathematics Institute (Tao Reference Tao2015). The rate of energy dissipation may not encapsulate all important features of the system even in the linear case. For example, the Coriolis force does not contribute to energy dissipation in general. Therefore, the recovery of the Stokes system, the Darcy equation and the BJSJ interface condition from the Helmholtz minimal dissipation principle is meaningful. In addition, the balance of the forces normal to the interface is nowhere used in the setup of the dissipation function. Hence it is a truly derived interface condition.

There are a variety of derivations of the Stokes–Darcy system, see for example the paper by Mikelic & Jäger (Reference Mikelic and Jäger2000). In fact, a two-phase version of the Stokes–Darcy system, the Cahn–Hilliard–Stokes–Darcy system, was derived via Onsager's extremum principle by Han, Sun & Wang (Reference Han, Sun and Wang2013) without much justification on the dissipation rate functions. See also the papers by Qian, Wang & Sheng (Reference Qian, Wang and Sheng2006) for a variational derivation of the moving contact line dynamics, and by Wang (Reference Wang2021) for a more general framework dubbed ‘generalized Onsager's extremum principle’. The formal derivation presented here can be made rigorous via saddle point theory. Details will be reported elsewhere.

2.1. The Stokes–Darcy system via Helmholtz's principle

We assume a rectangular domain $\varOmega = [0,L]\times [0,H]$ for simplicity. The domain $\varOmega$ is split into two disjoint subdomains $\varOmega _{c}$ and $\varOmega _{m}$ by a smooth curved interface $\varGamma$, with $\boldsymbol n$ denoting the unit normal of $\varGamma$ pointing from the conduit $\varOmega _{c}$ (the pure fluid region) to porous media $\varOmega _{m}$, and $\boldsymbol {\tau }$ denoting the unit tangent vector to $\varGamma$. See figure 1 for an illustration. The three-dimensional case can be formulated analogously.

Figure 1. Domian $\varOmega$.

2.1.1. The time-independent case

Helmholtz's minimal dissipation principle states that for a steady flow in a viscous liquid, with the speeds of flow on the boundaries of the fluid being given as steady and small, the Stokes flow is the one that minimizes the kinetic energy dissipation rate among all incompressible flow configurations (Batchelor Reference Batchelor2000). In other words, the fluid is smart.

Preliminary considerations. Because we are assuming small speeds, both velocities are taken to be incompressible. Hence the normal velocity across the interface $\varGamma$ between the conduit and the porous media must be continuous to ensure mass conservation, i.e.

(2.1)\begin{equation} \boldsymbol{u}_c\cdot{\boldsymbol n} = \boldsymbol{u}_m\cdot{\boldsymbol n}, \end{equation}

where $\boldsymbol {u}_c$ is the velocity in the conduit while $\boldsymbol {u}_m$ is the Darcy velocity in the porous media.

The rate of energy dissipation functions. To invoke Helmholtz's principle for the recovery of the coupled free flow and porous media flow systems, we identify the rate of kinetic energy dissipation. For the conduit part, the energy dissipation rate is $\int _{\varOmega _c} 2\eta _0 |\mathbb {D}(\boldsymbol {u}_c)|^2,$ where $\mathbb {D}(\boldsymbol {u}_c)=(\boldsymbol {\nabla } \boldsymbol {u}_c + \boldsymbol {\nabla } \boldsymbol {u}_c^T)/2$ is the rate of deformation tensor and $\eta _0$ is the (constant) viscosity of the fluid (Batchelor Reference Batchelor2000). For the porous media part, the energy dissipation rate is $\int _{\varOmega _m} ({\eta _0}/{\kappa }) |\boldsymbol {u}_m|^2,$ which can be inferred from the time-dependent Darcy equation (Le Bars & Worster Reference Le Bars and Worster2006; Nield & Bejan Reference Nield and Bejan2017). Here the permeability $\kappa$, a tensor in general, is taken to be homogeneous and isotropic, and hence can be represented as a positive constant for the sake of simplicity in our presentation. Our arguments remain valid in the general case as long as the heterogeneity in space is not too strong. Both the Stokes equation and the Darcy equation are used to motivate the rate of the kinetic energy dissipation function. However, their exact forms are not employed here. We reiterate that the rate of energy dissipation may not carry all the information of the underlying model.

Because there is an interface between the conduit and the porous media, additional dissipation on the boundary is possible. Indeed, despite the continuity of the normal velocities, the tangential velocities may contain a gap across the interface $\varGamma$. In their seminal work (Beavers & Joseph Reference Beavers and Joseph1967), Beavers and Joseph postulated that the viscous fluid exerts a force tangential to the interface (along the interface) that is proportional to the jump in tangential velocity $\boldsymbol {u}_c\cdot \boldsymbol {\tau } -\boldsymbol {u}_m\cdot \boldsymbol {\tau }$ to impede such a jump. A simple dimensional consideration suggests that this force takes the form $-({\alpha _0\eta _0}/{\sqrt {\kappa }})(\boldsymbol {u}_c\cdot \boldsymbol {\tau } -\boldsymbol {u}_m\cdot \boldsymbol {\tau })$, where $\alpha _0$ is a dimensionless number usually determined by experiment. Saffman (Reference Saffman1971) observed that the tangential component of the Darcy velocity at the interface $\varGamma$ is usually much smaller than its free flow counterpart, i.e. $|\boldsymbol {u}_c\cdot \boldsymbol {\tau }| \gg |\boldsymbol {u}_m\cdot \boldsymbol {\tau }|$, hence we may ignore $\boldsymbol {u}_m\cdot \boldsymbol {\tau }$ and approximate the force by $-({\alpha _0\eta _0}/{\sqrt {\kappa }})\boldsymbol {u}_c\cdot \boldsymbol {\tau }$. Therefore, the rate of work done on the interface to dissipate energy is approximately $\int _{\varGamma }({\alpha _0\eta _0}/{\sqrt {\kappa }})|\boldsymbol {u}_c\cdot \boldsymbol {\tau }|^2,$ where we have applied Saffman's simplification a second time.

In summary, the total kinetic energy dissipation rate takes the form:

(2.2)\begin{equation} F(\boldsymbol{u}_c,\boldsymbol{u}_m)=\int_{\varOmega_c} 2\eta_0 |\mathbb{D}(\boldsymbol{u}_c)|^2+ \int_{\varOmega_m} \frac{\eta_0}{\kappa} |\boldsymbol{u}_m|^2 +\int_{\varGamma} {\frac{\alpha_0\eta_0}{\sqrt{\kappa}}} |\boldsymbol{u}_{c}\cdot\tau|^2. \end{equation}

Boundary conditions. A variety of physically relevant boundary conditions are available. For simplicity, we impose a prescribed fluid velocity on $\varGamma _c$, the boundary of the conduit that is not shared with the porous media. The normal component is set to zero corresponding to no-forced filtration unless periodicity is assumed. And we impose a prescribed zero normal velocity on $\varGamma _m$, the boundary of the porous media not shared with the conduit unless suitable periodicity is stipulated.

The Lagrangian multipliers for the incompressibility constraints. To deal with the incompressibility constraints, we introduce two Lagrange multipliers $q_c, q_m$ for the incompressibility, and we may formulate Helmholtz's minimal dissipation principle for flows in pure fluid adjacent to porous media by minimizing the following functional:

(2.3)\begin{align} {\mathcal{F}}(\boldsymbol{u}_c, q_c, \boldsymbol{u}_m, q_m) &= \int_{\varOmega_c} 2\eta_0 |\mathbb{D}(\boldsymbol{u}_c)|^2- \int_{\varOmega_c} 2q_c \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_c +\int_{\varGamma} \frac{\alpha_0\eta_0}{\sqrt{\kappa}}|\boldsymbol{u}_c\cdot\boldsymbol{\tau}|^2 \nonumber\\ &\quad + \int_{\varOmega_m} \frac{\eta_0}{\kappa} |\boldsymbol{u}_m|^2 - \int_{\varOmega_m} 2q_m \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_m, \end{align}

with the boundary conditions specified above plus the mass conservation constraint $\boldsymbol u_m\cdot \boldsymbol n=\boldsymbol u_c\cdot \boldsymbol n$ on $\varGamma$. This formulation has the advantage that we do not need to deal with the incompressibility constraint explicitly, but through the Lagrange multipliers implicitly.

The derivation. Suppose that $(\boldsymbol {u}_c, {q}_c, {\boldsymbol {u}}_m, {q}_m)$ is a minimizer of the energy dissipation rate functional $\mathcal {F}$ with the constraint $\boldsymbol u_m\cdot \boldsymbol n=\boldsymbol u_c\cdot \boldsymbol n$ on $\varGamma$. Let $\varphi$ be a smooth function in $\varOmega _c$. Using the fact $(\boldsymbol {u}_c, q_c, \boldsymbol {u}_m, q_m)$ is a minimizer with the constraint $\boldsymbol u_m\cdot \boldsymbol n=\boldsymbol u_c\cdot \boldsymbol n$ on $\varGamma$, we deduce that the function

(2.4)\begin{equation} \varPhi_p(s)={\mathcal{F}}(\boldsymbol{u}_c, q_c+s\varphi, \boldsymbol{u}_m, q_m) \end{equation}

must attain its minimum at $s=0$. Hence, $\varPhi '_p(0)=0$, and we obtain:

(2.5)\begin{equation} - \int_{\varOmega_c} 2\varphi \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_c=0. \end{equation}

The above equation holds for any smooth function $\varphi$ on $\varOmega _{c}$. Hence we conclude that $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}_c=0$ in $\varOmega _{c}$. Similarly, $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}_m=0$ in $\varOmega _{m}$. Therefore, we recover the incompressibility.

Next, we consider perturbation in the velocity field. Let $\boldsymbol {v}_c$ be a smooth vector field on $\varOmega _c$ that vanishes on $\varGamma _c$, and $\boldsymbol {v}_m$ be a smooth vector field on $\varOmega _m$ satisfying $\boldsymbol {v}_m\cdot \boldsymbol n=0$ on $\varGamma _m$ and $\boldsymbol {v}_c\cdot \boldsymbol n=\boldsymbol {v}_m\cdot \boldsymbol n$ on the interface $\varGamma$. Then the function

(2.6)\begin{equation} \varPhi_v(s)={\mathcal{F}}(\boldsymbol{u}_c+s\boldsymbol{v}_c, q_c, \boldsymbol{u}_m+s\boldsymbol{v}_m, q_m) \end{equation}

attains its minimum at $s=0$ because $\boldsymbol {u}_c+s\boldsymbol {v}_c$ and $\boldsymbol {u}_m+s\boldsymbol {v}_m$ satisfy the specified boundary condition and the continuity of normal velocity constraint, and $(\boldsymbol {u}_c, {q}_c, {\boldsymbol {u}}_m, {q}_m)$ is a minimizer. Hence, $\varPhi '_v(0)=0$, which leads to

(2.7)\begin{align} & \int_{\varOmega_{c}}(4\eta_{0}\mathbb{D}(\boldsymbol{u}_{c}):\mathbb{D}(\boldsymbol{v}_c)-2q_c \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}_c})+ \int_{\varGamma}\frac{\alpha_{0}\eta_{0}}{\sqrt{\kappa}}2(\boldsymbol{u}_{c}\cdot{\boldsymbol\tau}) (\boldsymbol{v}_c\cdot{\boldsymbol\tau}) \nonumber\\ &\quad +\int_{\varOmega_{m}} \left(\frac{\eta_{0}}{\kappa}2\boldsymbol{u}_{m}\cdot\boldsymbol{v}_m-2q_m \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}_m\right)=0. \end{align}

Because $\mathbb {D}(\boldsymbol {u}_{c})$ is a symmetric matrix, it follows that $\mathbb {D}(\boldsymbol {u}_{c}):\mathbb {D}(\boldsymbol {v}_c)=\mathbb {D}(\boldsymbol {u}_{c}):\boldsymbol {\nabla }\boldsymbol {v}_c$. After performing integration by parts on the first and last terms on the left-hand side, we deduce

(2.8) \begin{align} 0 &= \int_{\varOmega_{c}}2({-}2\eta_{0}\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbb{D}(\boldsymbol{u}_{c})+\boldsymbol{\nabla} q_c)\cdot{\boldsymbol{v}_c} \nonumber\\ &\quad + \int_{\varGamma}\left(4\eta_{0}{\boldsymbol{v}_c}\cdot\mathbb{D}(\boldsymbol{u}_{c}){\boldsymbol n} +\frac{\alpha_{0}\eta_{0}}{\sqrt{\kappa}}2(\boldsymbol{u}_{c}\cdot{\boldsymbol\tau})({\boldsymbol{v}_c}\cdot{\boldsymbol\tau})\right) +\int_\varGamma 2(q_m \boldsymbol{v}_m\cdot\boldsymbol n - q_c \boldsymbol{v}_c\cdot\boldsymbol n) \nonumber\\ &\quad +\int_{\varOmega_{m}}2(\frac{\eta_{0}}{\kappa}\boldsymbol{u}_{m} + \boldsymbol{\nabla} q_m )\cdot{\boldsymbol{v}_m} \nonumber\\ &= \int_{\varOmega_{c}}2({-}2\eta_{0}\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbb{D}(\boldsymbol{u}_{c})+\boldsymbol{\nabla} q_c)\cdot{\boldsymbol{v}_c} + \int_{\varGamma}2(2\eta_{0}\boldsymbol n\cdot\mathbb{D}(\boldsymbol{u}_{c}){\boldsymbol n} -q_c+q_m)\boldsymbol{v}_c\cdot\boldsymbol n \nonumber\\ &\quad +\int_{\varGamma}2(2\eta_{0}{\boldsymbol\tau}\cdot\mathbb{D}(\boldsymbol{u}_{c}){\boldsymbol n} +\frac{\alpha_{0}\eta_{0}}{\sqrt{\kappa}}\boldsymbol{u}_{c}\cdot{\boldsymbol\tau})({\boldsymbol{v}_c}\cdot{\boldsymbol\tau}) +\int_{\varOmega_{m}}2\left(\frac{\eta_{0}}{\kappa}\boldsymbol{u}_{m} + \boldsymbol{\nabla} q_m \right)\cdot{\boldsymbol{v}_m}, \end{align}

where we have used the continuity of the normal velocity constraint $\boldsymbol {v}_c\cdot \boldsymbol n=\boldsymbol {v}_m\cdot \boldsymbol n$ on the interface $\varGamma$ in deducing the second term in the last equation.

Setting $\boldsymbol {v}_m=0$ and choosing $\boldsymbol {v}_c$ with $\boldsymbol {v}_c=0$ on $\varGamma$, we deduce

(2.9)\begin{equation} \int_{\varOmega_{c}}({-}2\eta_{0}\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbb{D}(\boldsymbol{u}_{c})+\boldsymbol{\nabla} q_c)\cdot{\boldsymbol{v}_c}=0, \end{equation}

which yields the Stokes equations

(2.10)\begin{equation} -2\eta_{0}\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbb{D}(\boldsymbol{u}_{c})+\boldsymbol{\nabla} q_c=0\quad\text{in}\ \varOmega_{c}, \end{equation}

with the Lagrangian multiplier $q_c$ for the incompressibility in the conduit serving the role of the pressure in the conduit.

Next, we set $\boldsymbol {v}_c=0$, and we deduce the Darcy equation

(2.11)\begin{equation} \frac{\eta_{0}}{\kappa}\boldsymbol{u}_{m}+\boldsymbol{\nabla} q_{m} =0, \end{equation}

with the Lagrangian multiplier $q_m$ for the incompressibility of the fluid in porous media serving the role of the pressure in porous media.

Hence we are left with

(2.12)\begin{align} 0 &= \int_{\varGamma}2(2\eta_{0}\boldsymbol n\cdot\mathbb{D}(\boldsymbol{u}_{c}){\boldsymbol n} -q_c+q_m)\boldsymbol{v}_c\cdot\boldsymbol n \nonumber\\ &\quad +\int_{\varGamma}2\left(2\eta_{0}{\boldsymbol\tau}\cdot\mathbb{D}(\boldsymbol{u}_{c}){\boldsymbol n}+ \frac{\alpha_{0}\eta_{0}}{\sqrt{\kappa}}\boldsymbol{u}_{c}\cdot{\boldsymbol\tau}\right)({\boldsymbol{v}_c}\cdot{\boldsymbol\tau}). \end{align}

Choosing $\boldsymbol {v}_c$ so that $\boldsymbol {v}_c\cdot \boldsymbol n=0$ on $\varGamma$, we deduce

(2.13)\begin{equation} \int_{\varGamma}\left(2\eta_{0}{\boldsymbol\tau}\cdot\mathbb{D}(\boldsymbol{u}_{c}){\boldsymbol n}+ \frac{\alpha_{0}\eta_{0}}{\sqrt{\kappa}}(\boldsymbol{u}_{c}\cdot{\boldsymbol\tau})\right){\boldsymbol{v}_c}\cdot{\boldsymbol\tau}=0. \end{equation}

Because ${\boldsymbol {v}_c}\cdot {\boldsymbol \tau }$ is arbitrary, from the equation above, we recover the BJSJ interface boundary condition:

(2.14)\begin{equation} {\boldsymbol\tau}\cdot\mathbb{T}(\boldsymbol{u}_{c}, q_c){\boldsymbol n}+ \frac{\alpha_{0}\eta_{0}}{\sqrt{\kappa}}\boldsymbol{u}_{c}\cdot{\boldsymbol\tau}=0. \end{equation}

On the other hand, choosing $\boldsymbol {v}_c$ so that $\boldsymbol {v}_c\cdot {\boldsymbol \tau }=0$ on $\varGamma$, we deduce

(2.15)\begin{equation} \int_{\varGamma}(2\eta_{0}\boldsymbol n\cdot\mathbb{D}(\boldsymbol{u}_{c}){\boldsymbol n} -q_c+q_m)\boldsymbol{v}_c\cdot\boldsymbol n=0. \end{equation}

This implies the balance of forces normal to the interface because $\boldsymbol {v}_c\cdot \boldsymbol n$ can be an arbitrary function,

(2.16)\begin{equation} \boldsymbol n\cdot\mathbb{T}(\boldsymbol{u}_c,q_c)\boldsymbol n ={-}q_m. \end{equation}

Summary. To summarize, we have derived the following Euler–Lagrange equations satisfied by the minimizer,

(2.17)$$\begin{gather} -2\eta_0\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbb{D}(\boldsymbol{u}_c) + \boldsymbol{\nabla} q_c = 0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_c=0,\quad \mbox{in}\ \varOmega_c, \end{gather}$$

(2.18)$$\begin{gather}\frac{\eta_0}{\kappa}\boldsymbol{u}_m +\boldsymbol{\nabla} q_m = 0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_m=0, \quad \mbox{in}\ \varOmega_m, \end{gather}$$

(2.19)$$\begin{gather}-{\boldsymbol\tau}\cdot \mathbb{T}(\boldsymbol u_c, q_c)\boldsymbol n = \frac{\alpha_0\eta_0}{\sqrt{\kappa}}{\boldsymbol\tau}\cdot\boldsymbol u_c, \quad \mbox{on}\ \varGamma, \end{gather}$$

(2.20)$$\begin{gather}-\boldsymbol n\cdot \mathbb{T}(\boldsymbol u_c, q_c)\boldsymbol n =q_m, \quad \mbox{on}\ \varGamma, \end{gather}$$

(2.21)$$\begin{gather}\boldsymbol u_c\cdot\boldsymbol n = \boldsymbol u_m\cdot\boldsymbol n, \quad \mbox{on}\ \varGamma. \end{gather}$$

Notice that as the last interfacial boundary condition, the continuity of normal velocity is imposed to ensure mass conservation. We observe that this Euler–Lagrange equation is exactly the classical Stokes–Darcy system with the BJSJ interface boundary condition (2.19) and the balance of the forces normal to the interface boundary condition (2.20), as proposed by Quarteroni and collaborators (Discacciati & Quarteroni Reference Discacciati and Quarteroni2009). The Lagrange multipliers for the incompressibility happen to be the pressures (hence the mathematical saying that ‘the pressure is the Lagrangian multiplier for the incompressibility’). Therefore, we will replace $q_c, q_m$ by $p_c, p_m$, respectively. Both the BJSJ interface boundary condition and the balance of the force normal to the interface $\varGamma$ appear in the process of this variational manipulation.

Remark 2.1 The presentation above indicates that the BJSJ interface boundary condition is fully consistent with the Helmholtz minimal dissipation principle. This is in accordance with the known fact that BJSJ leads to a well-posed Stokes–Darcy system with unique solution and continuous dependence on data. In contrast, there is no known direct consistency argument between the original BJ interface boundary condition and Helmholtz's minimal dissipation principle. There is no known well-posedness result for the steady Stokes–Darcy system with the original BJ interface boundary condition in general either. All these suggests the BJSJ interfacial boundary condition is more natural from the energetic consideration, and more convenient from the mathematics perspective.

2.1.2. Variational derivation of the case with a source

If there is a source term or external forcing $\boldsymbol {f}_c$ such as the gravitational force or a pressure gradient in the conduit, we could include an additional term related to the work done by the external forcing and arrive at the following generalized Helmholtz principle:

(2.22)\begin{align} {\mathcal{F}}_f(\boldsymbol{u}_c, p_c, \boldsymbol{u}_m, p_m) &= \int_{\varOmega_c} 2\eta_0 |\mathbb{D}(\boldsymbol{u}_c)|^2- \int_{\varOmega_c} 2p_c \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_c +\int_{\varGamma} \frac{\alpha_0\eta_0}{\sqrt{\kappa}}|\boldsymbol{u}_c\cdot\boldsymbol{\tau}|^2 \nonumber\\ &\quad + \int_{\varOmega_m} \frac{\eta_0}{\kappa} |\boldsymbol{u}_m|^2 - \int_{\varOmega_m} 2p_m \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_m -2\int_{\varOmega_c}\boldsymbol{f}_c\cdot\boldsymbol{u}_c, \end{align}

subject to the mass conservation constraint $\boldsymbol u_m\cdot \boldsymbol n=\boldsymbol u_c\cdot \boldsymbol n$ on $\varGamma$.

The same argument as in the previous subsection leads to the corresponding Euler–Lagrange equations for this case:

(2.23)$$\begin{gather} -2\eta_0\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbb{D}(\boldsymbol{u}_c) + \boldsymbol{\nabla} p_c = \boldsymbol{f}_c,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_c=0,\quad \mbox{in}\ \varOmega_c, \end{gather}$$

(2.24)$$\begin{gather}\frac{\eta_0}{\kappa}\boldsymbol{u}_m +\boldsymbol{\nabla} p_m = 0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_m=0, \quad \mbox{in}\ \varOmega_m, \end{gather}$$

together with the same interface coupling conditions (2.19)–(2.21). This is exactly the well-known time-independent (steady-state) Stokes–Darcy system with a source term in the conduit together with the BJSJ interface boundary condition and the balance of normal-force interface boundary condition (Discacciati et al. Reference Discacciati, Miglio and Quarteroni2002; Layton et al. Reference Layton, Schieweck and Yotov2002; Discacciati & Quarteroni Reference Discacciati and Quarteroni2009; Wilbrandt Reference Wilbrandt2019).

2.1.3. The time-dependent case

The time-dependent Stokes–Darcy system with BJSJ interface can be recovered analogously by using the gradient flow idea. The Darcy velocity $\boldsymbol {u}_m$ (or the pressure in porous media $p_m$) can be viewed as a function of the velocity $\boldsymbol {u}_c$ and the pressure $p_c$ in the conduit. Indeed, we could use (2.20) to solve Darcy's equation together with other appropriate boundary conditions once $\boldsymbol {u}_c$ and $p_c$ are given. In other words, the Darcy velocity is slaved by the Stokes velocity. Hence, we may formulate the objective functional as a functional of the variables in the conduit only and apply the gradient flow framework to come up with the following time-dependent Stokes–Darcy system:

(2.25)$$\begin{gather} \rho_0\frac{\partial\boldsymbol{u}_c}{\partial t}-\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbb{T}(\boldsymbol{u}_c, p_c) = \boldsymbol{f}_c,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_c=0,\quad \mbox{in}\ \varOmega_c, \end{gather}$$

(2.26)$$\begin{gather}\frac{\eta_0}{\kappa}\boldsymbol{u}_m +\boldsymbol{\nabla} p_m = 0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_m=0, \quad \mbox{in}\ \varOmega_m, \end{gather}$$

together with the same interface coupling conditions (2.19)–(2.21).

This is exactly the well-known time-dependent (unsteady-state) Stokes–Darcy system with a source term in the conduit together with the BJSJ interface boundary condition (Discacciati et al. Reference Discacciati, Miglio and Quarteroni2002; Discacciati & Quarteroni Reference Discacciati and Quarteroni2009; Wilbrandt Reference Wilbrandt2019).

Notice that the solution to this system is not unique. Indeed, if $(\boldsymbol {u}_c, p_c, \boldsymbol {u}_m, p_m)$ is a solution, $(\boldsymbol {u}_c, p_c+C, \boldsymbol {u}_m, p_m+C)$ is also a solution for any constant $C$. To ensure uniqueness of solution, we impose

(2.27)\begin{equation} \int_{\varOmega_m} p_m = 0. \end{equation}

This is equivalent to choosing a proper reference frame for the pressure.

2.2. Non-dimensional form

We now derive the non-dimensionalized form of the Stokes–Darcy system using the typical units associated with the free flow.

We use the typical length $L$, velocity $U_{0}$ and pressure difference $\delta P$ of the free flow in the conduit (pure fluid region) and we introduce the following scalings to non-dimensionalize the Stokes–Darcy system (2.25) and (2.26):

(2.28a–g)\begin{equation} \left.\begin{gathered} {\boldsymbol u}_{c} =\widetilde{\boldsymbol u}_{c}U_{0},\quad p_{c} =\tilde{p}_{c}\delta P,\quad {\boldsymbol x} =\widetilde{\boldsymbol x}L,\quad t =\tilde{t}\frac{L}{U_{0}},\\ {\boldsymbol f}_c =\widetilde{\boldsymbol f}F_{0},\quad {\boldsymbol u}_{m} =\widetilde{\boldsymbol u}_{m}U_{0},\quad p_{m} =\tilde{p}_{m}\delta P . \end{gathered}\right\} \end{equation}

This leads to

(2.29)\begin{equation} \left.\begin{gathered} \partial_{\tilde{t}}\widetilde{\boldsymbol u}_{c}-\frac{2}{Re}\widetilde{\boldsymbol{\nabla}}\boldsymbol{\cdot} \tilde{\mathbb{D}}(\widetilde{\boldsymbol u}_{c})+Eu\widetilde{\boldsymbol{\nabla}} \tilde{p}_{c} = \frac{Gr}{Re^{2}}\widetilde{\boldsymbol f},\quad \text{in}\ \varOmega_{c},\\ \widetilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\widetilde{\boldsymbol u}_{c} =0,\quad \text{in}\ \varOmega_{c}, \end{gathered}\right\} \end{equation}

and

(2.30)\begin{equation} \left.\begin{gathered} \frac{1}{Da\cdot Re}\widetilde{\boldsymbol u}_{m}={-}Eu\widetilde{\boldsymbol{\nabla}}\tilde{p}_{m},\quad \text{in}\ \varOmega_{m},\\ \widetilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\widetilde{\boldsymbol u}_{m}=0,\quad \text{in}\ \varOmega_{m}, \end{gathered}\right\} \end{equation}

where we have introduced the following four dimensionless parameters:

(2.31a–d)\begin{equation} Da=\frac{\kappa}{L^{2}}, \quad Re=\frac{\rho_{0} U_{0}L}{\eta_{0}}, \quad Gr=\frac{\rho_0 F_{0}L^{3}}{\eta_{0}^{2}}, \quad Eu=\frac{\delta P}{\rho_{0}U_{0}^{2}}, \end{equation}

which are the Darcy number, the Reynolds number, the generalized Grashof number and the Euler number, respectively. See the book by Foias et al. (Reference Foias, Manley, Rosa and Temam2001) for more on the generalized Grashof number.

On the interface $\varGamma$, we have

(2.32)\begin{equation} \left.\begin{gathered} \widetilde{\boldsymbol u}_{c}\cdot{\boldsymbol n} =\widetilde{\boldsymbol u}_{m}\cdot{\boldsymbol n}, \\ \boldsymbol{\tau}\cdot 2\tilde{\mathbb{D}}(\widetilde{\boldsymbol u}_{c}){\boldsymbol n} ={-} \frac{\alpha_{0}}{\sqrt{Da}}\widetilde{\boldsymbol u}_{c}\cdot\boldsymbol{\tau}, \\ {\boldsymbol n}\cdot\frac{2}{Re}\tilde{\mathbb{D}}(\widetilde{\boldsymbol u}_{c}) {\boldsymbol n}-Eu\,\tilde{p}_{c} ={-}Eu\,\tilde{p}_{m}. \end{gathered}\right\} \end{equation}

For applications such as flows in karst aquifer and hyporheic flows, the Reynolds number, Euler number and generalized Grashoff number are usually modest while the Darcy number is very small, of the order of $10^{-6}$ or smaller (Bear Reference Bear1988; Nield & Bejan Reference Nield and Bejan2017). For instance, for the laboratory set-up described in § 2.1 of the paper by Janssen et al. (Reference Janssen, Cardenas, Sawyer, Dammrich, Krietsch and de Beer2012), the permeability $\kappa = 1.5\times 10^{-11}\ \textrm {m}^2$, while the depth of the porous media is 9 cm, the width of the sand bed is 28 cm and the length of the sand bed is 1.5 m. Hence the Darcy number, which is defined as the ratio of the permeability to the typical length squared, is approximately $6.67\times 10^{-12}$. On the other hand, in the so-called ‘low-discharge’ case studied in their work, the mean horizontal free flow velocity is 0.07 m s$^{-1}$ and the Reynolds number is estimated at 1300, which is of the order of $10^3$. We also infer from figure 7 of the paper by Janssen et al. (Reference Janssen, Cardenas, Sawyer, Dammrich, Krietsch and de Beer2012) that the mean horizontal velocity in the porous media is approximately $2\ \textrm {cm}\ \textrm {hr}^{-1}\equiv 5.56 \times 10^{-6}$ m s$^{-1}$. Darcy's law implies that the horizontal pressure gradient in the porous media is approximately $3.71\times 10^{2}\ \textrm {kg}\ \textrm {m}^{-2}\ \textrm {s}^{-2}$. Hence the horizontal pressure drop in the porous media is approximately $5.56\times 10^2\ \textrm {kg}\ \textrm {m}^{-1}\ \textrm {s}^{-2}$. The horizontal pressure drop in the free flow is roughly the same by the approximate continuity of the pressure across the interface. Consequently, the Euler number is roughly 113.47. Therefore, the product of the Reynolds number and the Darcy number is of the order of $10^{-9}$ while the Euler number is of the order of $10^2$. Hence, it is reasonable to assume the smallness of $Re\, Da$ while holding the other parameters constant. The generalized Grashoff number can be taken to be zero because the only body force in this case is the gravitational force which bears no direct impact on the horizontal fluid motion, even in this turbulent case. For a laminar flow example, we follow the laboratory set-up presented in § 3 of the paper by Faulkner et al. (Reference Faulkner, Hu, Kish and Hua2009). The conductivity of the porous media is estimated as $6.19\times 10^{-2}\ \textrm {cm}\ \textrm {s}^{-1}$. This implies that the permeability is approximately $6.34\times 10^{-11}$. Using the length of the channel as the typical length, which is 47.8 cm, we deduce that the Darcy number is roughly $2.77\times 10^{-10}$. The measured outflow rate in the conduit/channel is approximately $2.70\ \textrm {ml}\ \textrm {s}^{-1} = 2.7\times 10^{-6}\ \textrm {m}^3\ \textrm {s}^{-1}$. Hence the horizontal velocity in the conduit is estimated to be $6.75\times 10^{-3}$ m s$^{-1}$. As a result, the corresponding Reynolds number is approximately$1.34\times 10^2$. Figure 6(a) in their work indicates that the head loss over a distance of 15 cm is approximately 0.08 cm. Therefore, the corresponding Euler number is approximately $1.72\times 10^2$. The generalized Grashoff number can be set to zero again because the only body force is the gravitational force which has no direct impact on horizontal flows. The product of the Reynolds and Darcy number is approximately $3.71\times 10^{-8}$ which is approximately five orders of magnitude smaller than the reciprocal of the Euler number. Therefore, it is reasonable for us to consider the small-Darcy-number regime while holding the other parameters fixed.

Denoting $\widetilde {\widetilde {p}}=Eu\,\tilde {p},\ {\widetilde {\widetilde {\boldsymbol f}}}= ({Gr}/{Re^2})\widetilde {\boldsymbol f}$, letting $\alpha =\alpha _{0}\sqrt {Re}$, $\varepsilon ^{2}=Da\,Re$, and dropping the tildes, one obtains the following non-dimensional Stokes–Darcy system:

(2.33)$$\begin{gather} \partial_{t}{\boldsymbol u}_{c}-\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbb{T}({\boldsymbol u}_{c},p_{c}) ={\boldsymbol f},\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol u}_{c} =0,\quad \text{in}\ \varOmega_{c}, \end{gather}$$

(2.34)$$\begin{gather}{\boldsymbol u}_{m}={-}\varepsilon^{2}\boldsymbol{\nabla} p_{m},\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol u}_{m}=0,\quad \text{in}\ \varOmega_{m}, \end{gather}$$

(2.35)$$\begin{gather}{\boldsymbol u}_{c}\cdot{\boldsymbol n} ={\boldsymbol u}_{m}\cdot{\boldsymbol n},\quad \text{on}\ \varGamma, \end{gather}$$

(2.36)$$\begin{gather}\boldsymbol{\tau}\cdot\mathbb{T}({\boldsymbol u}_{c},p_{c}){\boldsymbol n} ={-}\frac{\alpha}{\varepsilon}{\boldsymbol u}_{c}\cdot\boldsymbol{\tau}, \quad \text{on}\ \varGamma, \end{gather}$$

(2.37)$$\begin{gather}{\boldsymbol n}\cdot\mathbb{T}({\boldsymbol u}_{c},p_{c}){\boldsymbol n} ={-}p_{m}, \quad \text{on}\ \varGamma, \end{gather}$$

together with the initial condition:

(2.38)\begin{equation} {\boldsymbol u}_{c}|_{t=0} ={\boldsymbol u}_{0},\quad \text{in}\ \varOmega_{c}. \end{equation}

Here ${\boldsymbol u}_{c}$ is the non-dimensional velocity in the conduit, $p_c$ is the product of the Euler number and the non-dimensional pressure of free flow, $\mathbb {T}({\boldsymbol u}_{c},p_{c})=({2}/{Re})\mathbb {D}({\boldsymbol u}_{c})-p_{c}\mathbb {I}$, $\mathbb {D}({\boldsymbol u}_{c})={(\boldsymbol {\nabla }{\boldsymbol u}_{c}+\boldsymbol {\nabla }{\boldsymbol u}_{c}^{T})}/{2}$, ${\boldsymbol f}$ and ${\boldsymbol u}_{0}$ are the given non-dimensional external force and initial data, $\varepsilon ^{2}=Da\,Re$ is a small parameter, ${\boldsymbol u}_{m}$ is the non-dimensional Darcy velocity and $p_{m}$ is the product of the Euler number and the non-dimensional pressure of the fluid in the porous media.

The Darcy equation (2.34) can be reformulated in terms of the pressure (or the hydraulic head) as a Laplace equation after taking the divergence of (2.34)$_{1}$,

(2.39a,b)\begin{equation} {\rm \Delta} p_{m}=0 \quad \text{in}\ \varOmega_{m},\quad \int_{\varOmega_m} p_m=0. \end{equation}

We point out that different non-dimensionalizations are available, see for example the papers by Chen & Chen (Reference Chen and Chen1988) and Straughan (Reference Straughan2001) among others.

3. Asymptotic expansion and approximate solutions

With the small parameter $\varepsilon ^2=Da\,Re$ in hand, we employ expansion in $\varepsilon$ to obtain approximations with various orders of accuracy. The leading order is of particular importance because it coincides with the groundwater research community's heuristic semi-decoupled approach. The mathematical proof of the validity of the approximations will be furnished in Appendix A. A formal asymptotic expansion for the related Navier–Stokes–Darcy–Boussinesq system was carried out by one of the authors with collaborators using a different non-dimensionalization methodology (McCurdy, Moore & Wang Reference McCurdy, Moore and Wang2019).

We formally assume the following expansion:

(3.1)$$\begin{gather} {\boldsymbol u}_{j} ={\boldsymbol u}_{j}^{(0)}+\varepsilon {\boldsymbol u}_{j}^{(1)}+ \varepsilon^{2} {\boldsymbol u}_{j}^{(2)}+\cdots, \end{gather}$$

(3.2)$$\begin{gather}p_{j} =p_{j}^{(0)}+\varepsilon p_{j}^{(1)}+\varepsilon^{2} p_{j}^{(2)}+\cdots, \end{gather}$$

with the index $j=c$ or $m$.

Substituting this expression into the Stokes equations (2.33), the pressure formulation of the Darcy equation (2.39a,b) and the interface boundary conditions (2.35)–(2.37), we deduce the following.

1. (i) Leading order $O(1)$: Matching the $O(\varepsilon ^{0})$ terms yields

   (3.3)\begin{equation} \left.\begin{gathered} \partial_{t}{\boldsymbol u}_{c}^{(0)}-\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbb{T}({\boldsymbol u}_{c}^{(0)},p_{c}^{(0)})= {\boldsymbol f},\quad \text{in}\ \varOmega_{c}, \\ \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol u}_{c}^{(0)}=0,\quad \text{in}\ \varOmega_{c}, \\ {\boldsymbol u}_{c}^{(0)}=0,\quad \text{on}\ \varGamma, \\ {\boldsymbol u}_{c}^{(0)}|_{t=0}={\boldsymbol u}_{0},\quad \text{in}\ \varOmega_{c}. \end{gathered}\right\} \end{equation}
   and
   (3.4)\begin{equation} \left.\begin{gathered} {\rm \Delta} p_{m}^{(0)}=0,\quad \text{in}\ \varOmega_{m}, \\ p_{m}^{(0)}={-}{\boldsymbol n}\cdot\mathbb{T}({\boldsymbol u}_{c}^{(0)},p_{c}^{(0)}) {\boldsymbol n} =p_c^{(0)},\quad \text{on}\ \varGamma, \\ \int_{\varOmega_{m}}p_{m}^{(0)}=0, \\ {\boldsymbol u}_{m}^{(0)}=0,\quad \text{in}\ \varOmega_{m}, \end{gathered}\right\} \end{equation}
   where we have used the fact that $\boldsymbol {u}_c^{(0)}=0$ on $\varGamma$, i.e. (3.3)$_3$, to deduce the last equality in (3.4)$_2$. Indeed, thanks to the Galilean invariance of the system, a generic point on the interface can be taken as the origin and the tangent plane (line) is horizontal at this point without loss of generality. It is easy to check that $\boldsymbol n \cdot \mathbb {D}({\boldsymbol u}_c^{(0)})\boldsymbol n$ is equal to the partial derivative of the vertical velocity $w$ with respective to the vertical variable $z$. By incompressibility, it is equal to the negative of the sum of the partial derivative of $u$ with respect to $x$ and the partial derivative of $v$ with respect to $y$, where $u$ and $v$ are the first two (horizontal) components of the velocity while $x,y$ are the horizontal coordinates. The fact that $u$ and $v$ are zero on the interface and that the $x$-axis and $y$-axis are tangential to the interface at the origin implies that the $x$ and $y$ partial derivatives of $u$ and $v$ must be zero at the origin. Hence $\boldsymbol n \cdot \mathbb {D}({\boldsymbol u}_c^{(0)})\boldsymbol n=0$, which completes the proof of the continuity of pressure as the leading order behaviour at small Darcy numbers.

   Remark 3.1 An important observation is that the leading order dynamics is exactly the semi-decoupled approach advocated by some researchers from the groundwater studies community (Cardenas & Wilson Reference Cardenas and Wilson2007; Cardenas & Gooseff Reference Cardenas and Gooseff2008; Janssen et al. Reference Janssen, Cardenas, Sawyer, Dammrich, Krietsch and de Beer2012). The solution procedure is: (step 1) solve the free flow (Stokes) problem as if the porous media is not there at all, and hence no-slip on the interface $\varGamma$; (step 2) solve the Darcy equation with the pressure on the interface given by the pressure of the free flow (Stokes) flow in the conduit. The free flow pressure $p_c^{(0)}$ is completely determined when we enforce the average zero constraint (3.4)$_3$ for the pressure in porous media.

   Remark 3.2 Another interesting observation is that the continuity of the pressure (3.4)$_2$ as well as the continuity of tangential velocities (3.3)$_3$, (3.4)$_4$ are the leading order behaviours at small Darcy numbers of the coupled system recovered via the Helmholtz minimal dissipation principle.

2. (ii) First-order equations $O(\varepsilon )$: Matching the $O(\varepsilon )$ terms yields, for Stokes system,

   (3.5)\begin{equation} \left.\begin{gathered} \partial_{t}{\boldsymbol u}_{c}^{(1)}-\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbb{T}({\boldsymbol u}_{c}^{(1)},p_{c}^{(1)})=0,\quad \text{in}\ \varOmega_{c}, \\ \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol u}_{c}^{(1)}=0,\quad \text{in}\ \varOmega_{c}, \\ {\boldsymbol u}_{c}^{(1)}\cdot{\boldsymbol n}=0,\quad \text{on}\ \varGamma, \\ {\boldsymbol u}_{c}^{(1)}\cdot\boldsymbol{\tau}={-}\frac{1}{\alpha}\boldsymbol{\tau}\cdot \mathbb{T}({\boldsymbol u}_{c}^{(0)},p_{c}^{(0)}){\boldsymbol n},\quad \text{on}\ \varGamma, \\ {\boldsymbol u}_{c}^{(1)}|_{t=0}=0,\quad \text{in}\ \varOmega_{c}. \end{gathered}\right\} \end{equation}
   For the Darcy equation,
   (3.6)\begin{equation} \left.\begin{gathered} {\rm \Delta} p_{m}^{(1)}=0,\quad \text{in}\ \varOmega_{m}, \\ p_{m}^{(1)}={-}{\boldsymbol n}\cdot\mathbb{T}({\boldsymbol u}_{c}^{(1)},p_{c}^{(1)}){\boldsymbol n},\quad \text{on}\ \varGamma, \\ \int_{\varOmega_{m}}p_{m}^{(1)}=0, \\ \boldsymbol{u}_m^{(1)}=0. \end{gathered}\right\} \end{equation}

   Remark 3.3 Notice that the tangential component of the stress associated with the leading order Darcy flow is non-trivial in general, i.e. $\boldsymbol {\tau }\cdot \mathbb {T}({\boldsymbol u}_{c}^{(0)},p_{c}^{(0)}){\boldsymbol n}=\boldsymbol {\tau }\cdot \mathbb {D}({\boldsymbol u}_{c}^{(0)}){\boldsymbol n} \neq 0$. This implies, according to (3.5)$_4$, ${\boldsymbol u}_{c}^{(1)}\cdot \boldsymbol {\tau }\neq 0$. This further implies, together with (3.3)$_3$, that ${\boldsymbol u}_c\cdot \boldsymbol {\tau }=O(\sqrt {Da})$ on $\varGamma$. We also notice that $\boldsymbol {u}_m^{(1)}=0$ according to (3.6)$_4$. Hence the tangential velocities are discontinuous in general but the discontinuity is of the order of $\varepsilon =\sqrt {Da}$ although they are continuous at the leading order. Similarly, we usually have ${\boldsymbol n}\cdot \mathbb {D}({\boldsymbol u}_{c}^{(1)}){\boldsymbol n}\neq 0$. Consequently, $p_m^{(1)}\neq p_c^{(1)}$ in general. Hence, the pressure is usually discontinuous but the discontinuity is of the order of $\varepsilon \approx \sqrt {Da}$. Therefore, the seemingly contradicting interfacial conditions, such as the continuity and discontinuity of tangential velocity as outlined in the second paragraph of the paper by Le Bars & Worster (Reference Le Bars and Worster2006), the continuity of the pressure adopted by the groundwater studies community and the balance of normal forces to the interface, can be reconciled at the small-Darcy-number regime as an approximation of the ‘intrinsic’ interfacial boundary conditions recovered from the Helmholtz minimal dissipation principle.

   Same as above, we solve for $({\boldsymbol u}_{c}^{(1)},p_{c}^{(1)})$ from (3.5) first, then solve for $p_{m}^{(1)}$ from (3.6) with the help of the $({\boldsymbol u}_{c}^{(1)},p_{c}^{(1)})$ value on boundary $\varGamma$. The value of the pressure $p_{c}^{(1)}$ should be adjusted so that the mean zero constraint on the pressure in porous media is satisfied.

3. (iii) Matching the $O(\varepsilon ^{2})$ terms, we deduce from the Stokes equations,

   (3.7)\begin{equation} \left.\begin{gathered} \partial_{t}{\boldsymbol u}_{c}^{(2)}-\boldsymbol{\nabla}\boldsymbol{\cdot} \mathbb{T}({\boldsymbol u}_{c}^{(2)},p_{c}^{(2)})=0,\quad \text{in}\ \varOmega_{c}, \\ \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol u}_{c}^{(2)}=0,\quad \text{in}\ \varOmega_{c}, \\ {\boldsymbol u}_{c}^{(2)}\cdot{\boldsymbol n}={-}\boldsymbol{\nabla} p_{m}^{(0)}\cdot{\boldsymbol n},\quad \text{on}\ \varGamma, \\ {\boldsymbol u}_{c}^{(2)}\cdot\boldsymbol{\tau}={-}\frac{1}{\alpha}\boldsymbol{\tau}\cdot \mathbb{T}({\boldsymbol u}_{c}^{(1)},p_{c}^{(1)}){\boldsymbol n},\quad \text{on}\ \varGamma, \\ {\boldsymbol u}_{c}^{(2)}|_{t=0}=0,\quad \text{in}\ \varOmega_{c}. \end{gathered}\right\} \end{equation}
   and from the Darcy equation,
   (3.8)\begin{equation} \left.\begin{gathered} {\rm \Delta} p_{m}^{(2)}=0,\quad \text{in}\ \varOmega_{m}, \\ p_{m}^{(2)}={-}{\boldsymbol n}\cdot\mathbb{T}({\boldsymbol u}_{c}^{(2)},p_{c}^{(2)}){\boldsymbol n},\quad \text{on}\ \varGamma, \\ \int_{\varOmega_{m}}p_{m}^{(2)}=0,\\ {\boldsymbol u}_{m}^{(2)}={-}\boldsymbol{\nabla} p_{m}^{(0)},\quad \text{in}\ \varOmega_{m}. \end{gathered}\right\} \end{equation}

   Recall that $p_{m}^{(0)}$ is obtained from the leading order expansion (3.3). One can solve (3.7) first, and then (3.8) by adjusting the averaged value of the pressure $p_{c}^{(2)}$ so that the mean zero constraint on the pressure in porous media is fulfilled.

   Remark 3.4 We observe that this formal expansion implies that the Darcy velocity is of the order of $\varepsilon ^2=Da\,Re$ because $\boldsymbol {u}_m^{(0)}\equiv \boldsymbol {u}_m^{(1)}\equiv 0$ while $p_m^{(0)}$ is non-trivial in general. Hence the transportation in porous media occurs over a time scale proportional to $1/{Da\,Re}=1/\varepsilon ^2$, which is very long for small Darcy numbers and moderate Reynolds numbers as in the examples presented in § 2.1.1. Therefore, approximations, numerical schemes included, must be valid over this long-time scale to capture the transport behaviour.

   Remark 3.5 Another observation is that the asymptotic expansion is consistent with Saffman's heuristic simplification assumption (Saffman Reference Saffman1971). Indeed, because the Darcy velocity is of the order of $\varepsilon ^2=Da\,Re$ according to the previous remark, and because the tangential velocity of the free flow at the interface is of the order of $\varepsilon$ according to Remark 3.3, we see that $|\boldsymbol {u}_m\cdot \boldsymbol {\tau }| \ll |\boldsymbol {u}_c\cdot \boldsymbol {\tau }|$ on $\varGamma$ at the small-Darcy-number regime.

4. (iv) In general, for any $k\ge 2$ matching the $O(\varepsilon ^{k})$ terms yields, from the Stokes equations,

   (3.9)\begin{equation} \left.\begin{gathered} \partial_{t}{\boldsymbol u}_{c}^{(k)}-\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbb{T}({\boldsymbol u}_{c}^{(k)},p_{c}^{(k)})=0,\quad \text{in}\ \varOmega_{c}, \\ \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol u}_{c}^{(k)}=0,\quad \text{in}\ \varOmega_{c}, \\ {\boldsymbol u}_{c}^{(k)}\cdot{\boldsymbol n}={-}\boldsymbol{\nabla} p_{m}^{(k-2)}\cdot {\boldsymbol n},\quad \text{on}\ \varGamma, \\ {\boldsymbol u}_{c}^{(k)}\cdot\boldsymbol{\tau}={-}\frac{1}{\alpha}\boldsymbol{\tau}\cdot \mathbb{T}({\boldsymbol u}_{c}^{(k-1)},p_{c}^{(k-1)}){\boldsymbol n},\quad \text{on}\ \varGamma, \\ {\boldsymbol u}_{c}^{(k)}|_{t=0}=0,\quad \text{in}\ \varOmega_{c}, \end{gathered}\right\} \end{equation}
   and from the Darcy equation,
   (3.10)\begin{equation} \left.\begin{gathered} {\rm \Delta} p_{m}^{(k)}=0,\quad \text{in}\ \varOmega_{m}, \\ p_{m}^{(k)}={-}{\boldsymbol n}\cdot\mathbb{T}({\boldsymbol u}_{c}^{(k)},p_{c}^{(k)}){\boldsymbol n},\quad \text{on}\ \varGamma, \\ \int_{\varOmega_{m}}p_{m}^{(k)}=0, \\ \boldsymbol{u}_m^{(k)}={-}\boldsymbol{\nabla} p_m^{(k-2)}. \end{gathered}\right\} \end{equation}

We can now define approximate solutions of arbitrary order $k$:

(3.11)\begin{equation} \left.\begin{gathered} {\boldsymbol u}_{j}^{app,k} ={\boldsymbol u}_{j}^{(0)}+\varepsilon {\boldsymbol u}_{j}^{(1)}+ \varepsilon^{2} {\boldsymbol u}_{j}^{(2)}+\cdots +\varepsilon^k\boldsymbol{u}_j^{(k)},\\ p_{j}^{app,k} =p_{j}^{(0)}+\varepsilon p_{j}^{(1)}+\varepsilon^{2} p_{j}^{(2)}+\cdots + \varepsilon^k p_j^{(k)}, \end{gathered}\right\} \end{equation}

with $j=c$ or $m$.