Stokes–Darcy system, small-Darcy-number behaviour and related interfacial conditions

Abstract

We show that the Stokes–Darcy system, which governs flows through adjacent porous and pure-fluid domains in the two-domain approach without forced filtration, can be recovered from the Helmholtz minimal dissipation principle. While the continuity of normal velocity across the interface is imposed explicitly for mass conservation, only the Beavers–Joseph–Saffman–Jones (BJSJ) interface boundary condition is imposed implicitly, and the balance of the normal-force interface boundary condition appears naturally in the variational process. This set of interface boundary conditions is well-accepted in the mathematics community. We show that these interfacial boundary conditions, at the physically important small-Darcy-number regime, are consistent with continuity of pressure across the interface condition. The tangential velocity and pressure are discontinuous in general but the discontinuity is of the order of the square root of the Darcy number. Hence these interfacial conditions are all approximately consistent in the physically important small-Darcy-number regime. The leading order dynamics in the pure fluid zone is governed by the Stokes system with the no-slip no-penetration boundary condition on the interface between the free zone and porous media at a small Darcy number. The leading order non-trivial dynamics in porous media is governed by the Darcy equation with the pressure on the interface prescribed by the pressure of the leading order Stokes flow in the pure fluid zone. Such a semi-decoupled approach has long been used by the groundwater community. Our result is the first rigorous work quantifying the error of this intuitive approach and relating different interfacial conditions.

1. Introduction

The coupling and interaction of free flow with flows in porous media are ubiquitous. Well-known examples include flows in karst aquifers, hyporheic zones, air or oil filtration, proton-exchange membrane (PEM) fuel cells, as well as blood filtration in the human body. Hence it is important to study the coupled system.

When the flow speed is relatively low (small Reynolds number), the incompressible Stokes system is an excellent model for flow in the free zone (conduit). Time dependency could be included if the temporal evolution is important. Darcy's equation is a well-accepted model for flow in fluid-saturated porous media. Both systems are well-understood, see for example the works of Bear (1988), Tritton (1988), Batchelor (2000) and Temam (2001). The coupling of the two different systems at the interface is of great importance. We will focus on the less controversial case of no-forced filtration. However, there are still genuine disagreements among researchers on what are the ‘correct’ interface boundary conditions even in this case. Continuity of the normal velocity and the pressure are regarded as well-accepted by the fluid dynamics community according to Le Bars & Worster (2006). Nevertheless, the situation with other interfacial conditions is less clear. The second paragraph of the paper by Le Bars & Worster (2006) listed a few competing and sometimes contradicting interfacial conditions such as the Beavers–Joseph (BJ) or its simplification known as the Beavers–Joseph–Saffman–Jones (BJSJ) interface condition, the continuity or discontinuity of the tangential velocity and the continuity or discontinuity of the tangential shear stress among others. See also the works of Ochoa-Tapia & Whitaker (1995a), Ochoa-Tapia & Whitaker (1995b), Straughan (2008) and Zampogna & Bottaro (2016). In addition, Eggenweiler & Rybak (2020) and Rybak et al. (2020) have provided some evidence on the unsuitability of the BJ interface boundary for filtration problems (It seems that the evidence provided is consistent with the model error introduced by replacing the Stokes system on a perforated domain with the Darcy equation, see for example the work of Shen (2020). The error observed between the microscale model and the macroscale model might be explained in terms of the model error in the porous media without referring to the interfacial conditions.). On the other hand, there are three well-established interface boundary conditions within the mathematics community (Discacciati, Miglio & Quarteroni 2002; Layton, Schieweck & Yotov 2002) for the no-forced-filtration case: (i) the continuity of the normal velocity representing the conservation of mass; (ii) the linear balance of forces normal to the interface; and (iii) the BJSJ interface boundary condition. See the next section for more info. These interfacial conditions lead to the mathematical notion of ‘well-posedness’ of the coupled Stokes–Darcy system in the sense that for given reasonable data, there exists one and only one solution, and the solution depends on the data in a continuous fashion (small change to data leads to small change to the solution). The BJSJ interface boundary condition is a heuristic simplification of the original BJ interface boundary condition, which postulates that the tangential component of the normal stress in the free zone (conduit) at the interface is proportional to the jump in tangential velocity across the interface (Beavers & Joseph 1967). Beavers and Joseph's work is empirical, and the well-posedness, i.e. the existence and uniqueness of solution to the coupled Stokes–Darcy equations with the BJ condition has not been well-understood. (See the paper by Cao et al. (2010) for a related result.) It seems that two of the interfacial conditions are widely used by both the fluid dynamics and the mathematics communities: (i) the continuity of normal velocity across the interface, and (ii) the BJSJ interface condition. However, the two communities differ on the other interfacial condition. While the fluid dynamics and the groundwater studies communities adopt the continuity of pressure, the mathematics community embraces the balance of the forces normal to the interface, which implies the discontinuity of the pressure across the interface. Hence, there is a genuine need to clarify the relationship between different interfacial boundary conditions, even in this no-forced-filtration case.

Owing to the practical importance of the coupled problem, a lot of effort has been devoted to the development of accurate and efficient numerical methods for the coupled Stokes–Darcy system. See for example the papers by Discacciati et al. (2002), Layton et al. (2002), Discacciati & Quarteroni (2003), Chen et al. (2011), Chen et al. (2013), Chen et al. (2016) and Ervin et al. (2018), the survey paper by Discacciati & Quarteroni (2009), the monograph by Wilbrandt (2019), and the references therein for works from the mathematical side. A natural idea for efficiency is to design decoupled algorithms so that only one of the subproblems needs to be solved at each time step or iteration. The convergence and long-time stability of these decoupled schemes usually involve a stringent time step constraint for physically relevant small permeability (or hydraulic conductivity). On the other hand, some researchers in the groundwater studies community take a much simpler direct decomposition approach in the numerical simulation: one first solves the free flow problem (Stokes equation) with the no-slip boundary condition on the interface between the free zone (conduit, pure fluid region) and the porous media, followed by solving Darcy's equation with the pressure on the interface prescribed by the pressure of the solution of the Stokes problem on the interface, see for example the papers by Cardenas & Wilson (2007), Cardenas & Gooseff (2008) and Janssen et al. (2012). The groundwater research community's heuristic semi-decoupled approach is certainly very efficient because it involves the solution of two classical problems at once. A natural question is whether such a direct decoupling is valid. If not, what is the error.

To answer the two questions mentioned above, i.e. the validity of the groundwater research community's heuristic decoupled approach and the relationship between various interfacial boundary conditions, we first demonstrate that the interfacial conditions adopted by the mathematics community can be recovered from the classical Helmholtz minimal dissipation principle as long as appropriate dissipation functions are identified. We then non-dimensionalize the system and look at the physically important parameter regimes for possible simplification. The permeability of many common porous media is small. For example, the permeability of well-sorted sand and gravel is approximately $10^{-8}M^2$ or smaller (Bear 1988). Hence the Darcy number $Da$, which is defined as the ratio of the permeability to the typical length squared, is small for most reasonably sized domains. Therefore, it is of great physical relevance to study the small-Darcy-number asymptotic behaviour of the coupled Stokes–Darcy system.

Starting from the system recovered from the Helmholtz minimal dissipation principle, we discover that at the small-Darcy-number regime the leading order dynamics is exactly the semi-decoupled dynamics proposed by the groundwater studies community. We also discover the following relationships between different interfacial boundary conditions at the small-Darcy-number regime: (i) the continuity of the pressure interface condition adopted by the fluid dynamics and groundwater studies communities in the no-forced-filtration case is the leading order behaviour of the balance of forces normal to the interface condition adopted by the mathematics community, and recovered via the minimal dissipation principle, see Remark 3.2; (ii) the continuity of tangential velocity across the interface is the leading order behaviour of the coupled system while the tangential velocity is discontinuous across the interface in general, but the discontinuity is of the order of the square root of the Darcy number, see Remark 3.3; (iii) Saffman's approximation is valid in the sense that the tangential velocity of the Darcy flow at the interface is of the order of the Darcy number, while the tangential velocity of the free flow is of the order of the square root of the Darcy number, see Remark 3.5.

The rest of the manuscript is organized as follows. In § 2, we present a derivation of the Stokes–Darcy system together with the ‘natural’ interface boundary conditions based on the Helmholtz minimal dissipation principle for the no-forced-filtration case. We then introduce the non-dimensional governing equations together with the boundary and interface conditions. In § 3, we construct approximate solutions via asymptotic expansion in the small parameter Darcy number. In particular, we show that the Stokes equation with the no-slip boundary condition on the interface is the leading order dynamics for the pure fluid region, and Darcy's equation with the pressure on the interface prescribed by the pressure of the Stokes flow governs the leading order behaviour in porous media. The consistency of various interfacial boundary conditions is discussed. A rigorous mathematical proof of the formal asymptotic behaviour on a time scale that is relevant to transport in porous media is included in Appendix A for interested readers.

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 2000) 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 2015). 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 (2000). 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 (2013) without much justification on the dissipation rate functions. See also the papers by Qian, Wang & Sheng (2006) for a variational derivation of the moving contact line dynamics, and by Wang (2021) 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 2000). 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 2000). 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 2006; Nield & Bejan 2017). 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 1967), 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 (1971) 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 2009). 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. 2002; Layton et al. 2002; Discacciati & Quarteroni 2009; Wilbrandt 2019).

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. 2002; Discacciati & Quarteroni 2009; Wilbrandt 2019).

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. (2001) 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 1988; Nield & Bejan 2017). For instance, for the laboratory set-up described in § 2.1 of the paper by Janssen et al. (2012), 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. (2012) 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. (2009). 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 (1988) and Straughan (2001) 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 2019).

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 2007; Cardenas & Gooseff 2008; Janssen et al. 2012). 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 (2006), 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 1971). 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$.

Remark 3.6 Heuristically, the error in approximating $\boldsymbol {u}$ by $\boldsymbol {u}^{app,k}$ is of the order of $\varepsilon ^{k+1}$. The semi-decoupled approach proposed by the groundwater research community corresponds to $k=0$. Hence the error is of the order of $\varepsilon =\sqrt {Da}\sqrt {Re}$ formally. Higher order approximation with smaller error can be achieved by including more terms in the expansion at the expense of solving a few more Stokes equations and Darcy equations. Alternatively, iterations can be used to reduce the error as well (Li et al. 2020).

Remark 3.7 The terms in the expansion could be made as smooth as we need as long as the initial data, the external forcing term and the interface are smooth enough and certain compatibility conditions are satisfied (Temam 1982). The terms would remain bounded in time provided that the external forcing term remains bounded in time in an appropriate fashion. The smoothness of the data, their compatibility and boundedness in time have been assumed throughout this manuscript.

4. Conclusion

We have recovered the Stokes–Darcy system via Helmholtz's minimal dissipation principle with appropriate dissipation functions for the no-forced-filtration case. The Stokes equation, the Darcy equation and the Beavers–Joseph interface boundary condition are used in motivating the dissipation functions, but not elsewhere in the derivation process. It turns out that the balance of the forces normal to the interface boundary condition emerge as ‘intrinsic’ in the sense that they appear naturally out of the minimal dissipation set-up. We also recover the celebrated Beavers–Joseph–Saffman(–Jones) interface boundary condition. The system is non-dimensionalized using the typical scales associated with the conduit (free flow zone, pure fluid zone). We performed an asymptotic expansion of the non-dimensional system with respect to the physically small parameter Darcy number. The leading order expansion recovers the solution offered by the groundwater studies community, i.e. solving the Stokes problem with no-slip boundary condition on the interface, and then the Darcy equation with the pressure on the interface prescribed by the pressure of the flow in the conduit (free flow). We also discover that the continuity of pressure and the continuity of the tangential velocities interfacial conditions are all consistent with the leading order behaviour of the Stokes–Darcy system at small Darcy numbers. The leading order Darcy velocity is of the order of the Darcy number. Hence transport in porous media occurs on a time scale of the order of the reciprocal of the Darcy number. Saffman's assumption is also recovered in the formal asymptotic expansion. A rigorous mathematical proof of the formal expansion valid over the physically important transportation time scale is furnished in Appendix A.

There are several directions that may merit future attention. The first is the nonlinear case with the Stokes equation replaced by the Navier–Stokes equation. A formal asymptotic expansion for the related Navier–Stokes–Darcy–Boussinesq system was carried out by one of the authors with collaborators albeit using a different non-dimensionalization methodology in the paper by McCurdy et al. (2019). See also the papers of Chen & Chen (1988) and Straughan (2001). However, no rigorous result is available so far. The second is the case with the original interface boundary condition proposed byBeavers & Joseph (1967) instead of the simplified Beavers–Joseph–Saffman–Jones interface boundary condition. The difficulty is partly associated with the well-posedness issue of the coupled system with the Beavers–Joseph interface boundary condition. A careful analysis should be able to establish the BJSJ as a good approximation of the BJ interface boundary condition. This would provide additional rigorous mathematical foundation for Saffman's argument (Saffman 1971). In the case of forced filtration, the BJ interface boundary condition may not be appropriate and hence the investigation of a suitable interface boundary condition is even more interesting (Eggenweiler & Rybak 2020; Rybak et al. 2020). The third is the development of accurate and efficient numerical methods that are valid over a long time (with an error estimate that is uniform for the time interval proportional to the reciprocal of the Darcy number for instance). The natural decoupling furnished by the asymptotic expansion should aid with the development of decoupled schemes that use the solver for the subproblems only and avoid the stringent time step restriction associated with those long-time accurate numerical schemes developed earlier, see for example the paper by Chen et al. (2016).