# Mathematical study of air flow models around wooded areas

Laurent Chupin, Marguerite Gisclon

<figure><img src="image_1.png" /></figure>

## To cite this version:

Laurent Chupin, Marguerite Gisclon. Mathematical study of air flow models around wooded areas. 2025. hal-04913287

HAL Id: hal-04913287

https://hal.science/hal-04913287v1

Preprint submitted on 27 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mathematical study of air flow models around wooded areas

Laurent Chupin,\* Marguerite Gisclon†

January 21, 2025

## Abstract

The air flow in a natural environment is generally disturbed by natural elements that are not always taken into account in mathematical models. In this article, we aim to propose a model that correctly accounts for the effect of a wooded area (e.g. trees, hedges or a forest) on the air flow. If the air flow in a "free" environment is well described by the Navier-Stokes equations, our aim will be to propose choices for the air flow in a wooded area. Typically, we have relied on models based on porous media, such as Darcy's or Brinkman's models. The key point is then to understand how to model the interface between these two domains, the "free" air and the "wooded" zone. We will propose several options, with the aim of ensuring that the links between the various choices are continuous and relatively easy to implement. In fact, the results we propose are both a mathematical analysis of the suggested models (existence of solution, continuity with respect to parameters), and results of numerical simulations highlighting the differences between the models.

## 1 Introduction

The aim of this article is to propose and implement a model that will enable us to numerically experiment with different planting scenarios in an urban environment, in order to assess their impact in terms of micro-climate improvement. The approach, based on mathematical modelling, is designed to provide versatility for numerical implementation. The objective is to describe the airflow through an urban area including wooded areas and settlements, combining porous and fluid regions, i.e. two models with interface conditions. As flow models are fairly classical (Stokes or Navier-Stokes in air, and Darcy in a porous medium, for example), the choice of a model at the interface is crucial. It is on this last point that we propose to establish a discussion and put forward the mathematical and numerical analysis of several models.

The rest of this introduction will be devoted to a review of the state of the art in coupling air flows and flows in porous media. The introduction concludes with a presentation of the geometry in which we have been working. Section 2 presents the notations and the Stokes model in the fluid part of the flow, while Section 3 focuses on the different models that can be proposed in the porous medium, taking interface choices into account. Three different models are proposed. These three models are mathematically analyzed in Section 4. We show that they are all well-posed (more precisely, that solutions exist and are unique). It is also shown that there is continuity between these models, and that it is possible to move from one model to another by varying the parameters. Section 5 proposes an extension of the model in air, taking non-linear terms into account (more precisely, using the Navier-Stokes equations instead of the Stokes equation to describe air flow). Finally, these different models and different geometries are highlighted in Section 6, where numerical simulations of the models presented are carried out. They show the influence of model and parameter choices on flow in variable geometries.

### 1.1 State of the art

There exist a number of works on coupling conditions at the sharp fluid-porous interface for the Stokes-Darcy models that are based on the Beavers-Joseph condition. The coupling is commonly modelled by the interface condition postulated by Beavers and Joseph in [9] or by its simplification introduced by Saffman in [27] and called the Beavers-Joseph interface condition. Many of them use a wide range of techniques: asymptotic modelling in [5] and [6], homogenisation methods in [18] and [28]. For example, in [28], the authors consider steady-state fluid flow in the coupled domain consisting of a free-flow region and a porous medium and reformulate the generalized interface condition on the tangential component of velocity for arbitrary flows in Stokes/Darcy systems, so that it has the same analytical form as the Beavers-Joseph condition. The porous layer contains periodically distributed solid inclusions that enable the permeability and effective boundary condition to be calculated. Other techniques are penalization theory in [20], regularization in [12], boundary layer theory in [8], [14] and [18], or upscaling method as in [7]. In [7], the authors derive an original nonlinear multi-dimensional model for the inertial fluid flow through a fluid-porous interface by asymptotic theory for arbitrary flow directions. On the bottom of the transition zone, they consider the classical set of interface conditions based on the Beavers-Joseph-Saffman approach (see [16]).

<footnote>\*Laboratoire de Mathématiques Blaise Pascal, Université Clermont Auvergne</footnote>

<footnote>†Laboratoire de Mathématiques, Université Savoie Mont Blanc</footnote>
They prove that the both coupled problems are well posed.

From theoretical point of view, many authors have studied the coupling of the Stokes and Darcy systems to model the filtration of incompressible fluids through porous media, and we can only list very few of them. Here are some non-exhaustive works who analyze different global coupled model in order to prove its well-posedness. P. Angot and his collaborators have a lot of results about this subject. In [4], the authors establish rigorous estimates of the error induced by a $L^2$ penalization inducing in a Darcy equation in the solid body and by a $H^1$ penalization inducing in a Brinkman equation in the body. The solvability is proved in [2] with no restriction on the size of the slip coefficient. Next, in [5], the authors consider the coupling of a Brinkman model and Stokes equations with jump embedded transmission conditions by assuming that the viscosity in the porous region is very small.

More recently, other studies have focused on these issues. In [8], the authors use the porosity of the media as optimization parameter hence to minimize some cost function by finding the location of a porous media and show that the weak-limit is a fluid-porous media interface problem given by a penalization model for Navier-Stokes-Darcy problem. In [21], the authors consider a mixed model for coupling a fluid flow and a porous media flow in a bounded smooth domain. Their interface conditions are the mass conservation, the balance of normal force and the Beavers-Joseph interface condition. This last condition means the tangential components of the normal stress force are proportional to the difference of the tangential components of the fluid flow and the porous media flow velocities. They get the uniqueness of the possible solution to the steady-state coupled Stokes/Darcy model and they show the existence of solutions by using Galerkin method. In [20], the authors incorporate rigid moving obstacles in a fluid dynamics context using concepts from porous media theory. Based on the Navier-Stokes-Brinkman equations which augments the Navier-Stokes equation with a Darcy drag term their method, using a Navier-Stokes-Brinkman penalization approach, represents solid obstacles as time-varying regions containing a porous medium of vanishing permeability.

In this article, the fluid flow is not necessarily parallel to the fluid-porous interface as in the work of P. Angot et al. ([6] and [3]) where it is relatively easy to consider a transition region and thus propose Stokes-Brinkman-Darcy models with appropriate coupling conditions at the top and bottom. In our work, this requires taking into account a Beavers-Joseph coupling condition typically used for the Stokes-Darcy problem in its entirety.

## 1.2 Geometric configuration

From a mathematical point of view, we are interested in an air flow in a natural environment represented by an open boundary $\Omega$ of $\mathbb{R}^2$, whose boundary is denoted $\Gamma_f = \partial\Omega$. Within this domain, an open $\Omega_p$ represents the porous parts (woods, trees, forests...). We note $\Gamma_p$ the boundary of $\Omega_p$. The complementary domain $\Omega_f = \Omega \setminus \bar{\Omega}_p$ thus corresponds to a fluid zone in which air can circulate freely. In terms of regularity, we will assume that these domains are Lipschitz.

Designate by $\mathbf{n}$ the unit normal vector at the $\Gamma_p$ interface, emerging from the $\Omega_p$ domain. The vector $\mathbf{t}$ corresponds to the unit tangent vector at the interface $\Gamma_p$, so that the family $(\mathbf{n}, \mathbf{t})$ forms a direct basis of the plane. Figure 1 summarizes these notations and the proposed geometric configuration.

<figure><img src="image_3.png" /><figcaption>Figure 1: Geometry and notations</figcaption></figure>

**Remark 1** The case of an open $\mathbb{R}^3$ is similar, but more technical since, in particular, we have to define a tangent plane to a surface rather than a tangent. The mathematical analysis in this three-dimensional case would be unchanged. The only dimension-dependent ingredients are the Sobolev's injections, but those used in this article and therefore in the two-dimensional case are still valid in the three-dimensional case.
## 2 Flow in the fluid domain

**Equation inside the flow** The flow outside the wooded areas is a simple air flow assumed to be non-turbulent. We then assume that the velocity $\mathbf{u}_f$ (m.s$^{-1}$) and the pressure $p_f$ (Pa) in the domain $\Omega_f$ satisfy the Stokes's system

$$ \left\{ \begin{array}{l} -\eta \Delta \mathbf{u}_f + \nabla p_f = \mathbf{F} \quad \text{in } \Omega_f \\ \text{div } \mathbf{u}_f = 0 \quad \text{in } \Omega_f \end{array} \right. \qquad (1) $$

$$ \text{div } \mathbf{u}_f = 0 \quad \text{in } \Omega_f \qquad (2) $$

The coefficient $\eta > 0$ corresponds to the viscosity of the air and for numerical applications we will use the value $\eta = 2 \times 10^{-5}$ Pa.s. The source term $\mathbf{F}$ represents any external forces taken into account.

**Remark 2** We chose to use the Stokes's model which takes into account neither temporal nor inertial effects. Indeed, we wish to capture phenomena on time scales long enough to neglect the term $\partial_t \mathbf{u}_f$, and in a laminar regime that does not write the term $\mathbf{u}_f \cdot \nabla \mathbf{u}_f$. In practice, it is entirely possible to take this nonlinear term into account in the mathematical analysis. For the sake of clarity, we did not write it afterwards but Section 5 is entirely dedicated to adding this term.

For the following, it's important to notice that we can write the Stokes's equation (1) using the stress tensor $\sigma_f = 2\eta \mathbf{D} \mathbf{u}_f - p_f \mathbf{I} \mathbf{d}$, the notation $\mathbf{D} \mathbf{u} = \frac{1}{2}(\nabla \mathbf{u} + (\nabla \mathbf{u})^\top)$ standing for the symmetric part of the velocity gradient. In terms of stress, the equation (1) writes

$$ - \text{div} \sigma_f = \mathbf{F} \quad \text{in } \Omega_f \qquad (3) $$

**Condition on the outer boundary** To correctly describe a situation based on Stokes equations, we need boundary conditions. We will impose homogeneous Dirichlet conditions on the "outer" boundary:

$$ \mathbf{u}_f = \mathbf{0} \quad \text{on } \Gamma_f \qquad (4) $$

**Remark 3** Notice that it is possible to consider non-homogeneous boundary conditions

$$ \mathbf{u}_f = \mathbf{u}_{\text{bound}} \quad \text{on } \Gamma_f. $$

The only restriction is linked to the incompressibility condition (2): we must have $\int_{\Gamma} \mathbf{u}_{\text{bound}} \cdot \mathbf{n} = 0$.

From a mathematical point of view, in the case of a non-homogeneous Dirichlet condition, it is not difficult to return to the homogeneous case, even if it requires modifying the source term $\mathbf{F}$. To do this, we extend $\mathbf{u}_{\text{bound}}$ on $\Omega$ into a vector field $\tilde{\mathbf{u}}_{\text{bound}}$ with the regularity $H^1(\Omega)$, free-divergence and zero on $\Omega_p$ (this construction is possible, via Bogovskii's operator given in [10], see also [11]). The difference $\mathbf{u}_f - \tilde{\mathbf{u}}_{\text{bound}}$ vanishes on $\Gamma_f$ and satisfies the Stokes's problem (1)-(2) by replacing $\mathbf{F}$ by $\mathbf{F} + \eta \Delta \tilde{\mathbf{u}}_{\text{bound}}$.

**Condition on the inner boundary** The conditions on the inner boundary $\Gamma_p$, which corresponds to the interface between the air flow and the wooded areas, depend on the flow model chosen for the porous medium $\Omega_p$. Indeed, several models linking velocity and pressure in a porous medium are possible. We can use the Darcy's model or the Brinkman's model. In both cases, we need to specify the interface conditions on the boundary $\Gamma_p$.

In the next section we present different possible couplings.

# 3 Flow in the porous medium

## 3.1 Stokes-Darcy model

The "simplest" model to describe a flow in a porous medium $\Omega_p$ is the Darcy's model. We note $\mathbf{u}_p$ the velocity and $p_p$ the pressure in $\Omega_p$, then the Darcy's model writes

$$ \left\{ \begin{array}{l} \kappa \mathbf{u}_p + \nabla p_p = \mathbf{0} \quad \text{in } \Omega_p, \\ \text{div } \mathbf{u}_p = 0 \quad \text{in } \Omega_p. \end{array} \right. \qquad (5) $$

$$ \text{div } \mathbf{u}_p = 0 \quad \text{in } \Omega_p. \qquad (6) $$

The coefficient $\kappa > 0$ corresponds to the ratio between the viscosity $\eta$ (Pa.s) and the permeability $k$ (m$^2$) which measures the ability of the medium to allow the fluid to pass through it.

In the case of modeling the foliage of a tree or a wooded area that the domain would represent the domain $\Omega_p$, the Kozeny-Carman's law allows us to have an order of magnitude of the permeability and therefore of the coefficient $\kappa$. This law writes

$$ k = \frac{d^2 \phi^3}{150(1 - \phi)^2} \qquad (7) $$

where $\phi$ represents the porosity of the foliage and $d$ the equivalent diameter of an average leaf/branch. As an example (other examples will be given in the section 6 devoted to the numerical simulations, see for instance Table 1, page 12), if we choose $\phi = 0.9$ and $d = 0.05$ m we obtain

$$ k \approx 1.21 \times 10^{-3} \text{ m}^2 $$
then (we recall that $\eta = 2 \times 10^{-5}$ Pa.s)

$$ \kappa \approx 1.64 \times 10^{-2} \text{Pa.s.m}^{-2}. $$

To couple this model (5)-(6) with the Stokes's model (1)-(2) we impose conditions on the interface $\Gamma_p$. We use the Beavers-Joseph-Saffman's conditions, see the works of Beavers and Joseph [9] completed by those of Saffman [27] and reconsidered in the fundamental paper [19]. These are relationships coupling the velocity $\mathbf{u}_f$ and $\mathbf{u}_p$ but also $\boldsymbol{\sigma}_f$ and $\boldsymbol{\sigma}_p = -p_p \mathbf{Id}$.

To understand the interface conditions, we use the arguments of [22]

1. Mass conservation across $\Gamma_p$ is expressed by $(\mathbf{u}_f - \mathbf{u}_p) \cdot \mathbf{n} = 0$.

2. Continuity of forces gives $((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{n} = 0$.

3. The last boundary condition in best agreement with experimental evidence evolved from the work of Beavers and Joseph [9] and states that "(shear stress along $\Gamma_p$) is proportional to (slip velocity along $\Gamma_p$)". Mathematically, this can be represented by $((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{t} = \beta(\mathbf{u}_f - \mathbf{u}_p) \cdot \mathbf{t}$ where $\beta > 0$.

However, it has been observed that the term on the right-hand side $\mathbf{u}_p \cdot \mathbf{t}$ is much smaller than the other terms. The most accepted interface condition was derived by Saffman [27] using a statistical approach and the Brinkman approximation and also by Jäger and Mikelić [23] who gave a mathematical justification. Indeed, the reader can refer to the work of Jäger and Mikelić for the derivation by homogenization of the Beavers-Joseph-Saffman interface condition. This condition, which drops this term, is now known as the Beavers-Joseph-Saffman law and is thus given by $((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{t} = \beta \mathbf{u}_f \cdot \mathbf{t}$.

With these boundary conditions, the model writes

$$ \begin{array}{ll} \text{(SD)} & \left\{ \begin{array}{ll} -\eta \Delta \mathbf{u}_f + \nabla p_f = \mathbf{F} & \text{in } \Omega_f \\ \text{div } \mathbf{u}_f = 0 & \text{in } \Omega_f \\ \kappa \mathbf{u}_p + \nabla p_p = \mathbf{0} & \text{in } \Omega_p \\ \text{div } \mathbf{u}_p = 0 & \text{in } \Omega_p \\ (\mathbf{u}_f - \mathbf{u}_p) \cdot \mathbf{n} = 0 & \text{on } \Gamma_p \\ ((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{n} = 0 & \text{on } \Gamma_p \\ ((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{t} = \beta \mathbf{u}_f \cdot \mathbf{t} & \text{on } \Gamma_p \\ \mathbf{u}_f = \mathbf{0} & \text{on } \Gamma_f \end{array} \right. \end{array} \quad \begin{array}{l} (8) \\ (9) \\ (10) \\ (11) \\ (12) \\ (13) \\ (14) \\ (15) \end{array} $$

From a numerical value point of view the coefficient $\beta$ was obtained experimentally in [9]

$$ \beta = \frac{\eta \alpha}{\sqrt{k}} $$

where $\alpha$ is a dimensionless number which generally takes values of the order of 1 (see [9] or [24]). With the numerical values previously used we obtain $\beta \approx 5.77 \times 10^{-4}$ Pa.s.m$^{-1}$.

## 3.2 Stokes-Brinkman's model

It is also possible to use Brinkman's model to describe porous flow in $\Omega_p$. This model consists of including a diffusive part in the stress of the porous medium. More precisely Equation (10) is replaced by

$$ -\gamma \Delta \mathbf{u}_p + \kappa \mathbf{u}_p + \nabla p_p = \mathbf{0} \quad \text{in } \Omega_p \quad (16) $$

where $\gamma > 0$ quantifies diffusion and is homogeneous to a viscosity (Pa.s). This formulation is suitable for describing flow in a medium with a high porosity [17, 25, 26]. In the case of modelling a wooded region, the porosity seems very high (of the order of 0.9) which may justify this approach. However, obtaining a value for the effective viscosity $\gamma$ is tricky [25] and many authors choose $\gamma = \eta$, see for example [1, 26].

From a theoretical point of view, the additional diffusive term implies that the conditions (12)-(14) are insufficient to ensure the existence and uniqueness of a solution. We will then add the following condition

$$ (\mathbf{D}\mathbf{u}_p \cdot \mathbf{n}) \cdot \mathbf{t} = 0 \quad \text{on } \Gamma_p. \quad (17) $$

Making appear the stress $\boldsymbol{\sigma}_p = 2\gamma\mathbf{D}\mathbf{u}_p - p_p \mathbf{Id}$ this condition is equivalent to

$$ (\boldsymbol{\sigma}_p \cdot \mathbf{n}) \cdot \mathbf{t} = 0 \quad \text{on } \Gamma_p. \quad (18) $$
Finally, Stokes-Brinkman’s model writes

$$
\begin{equation}
\begin{aligned}
\text{(SB)} \quad \left\{
\begin{array}{ll}
-\eta \Delta \mathbf{u}_f + \nabla p_f = \mathbf{F} & \text{in } \Omega_f \\
\operatorname{div} \mathbf{u}_f = 0 & \text{in } \Omega_f \\
-\gamma \Delta \mathbf{u}_p + \kappa \mathbf{u}_p + \nabla p_p = \mathbf{0} & \text{in } \Omega_p \\
\operatorname{div} \mathbf{u}_p = 0 & \text{in } \Omega_p \\
(\mathbf{u}_f - \mathbf{u}_p) \cdot \mathbf{n} = 0 & \text{on } \Gamma_p \\
((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{n} = 0 & \text{on } \Gamma_p \\
(\boldsymbol{\sigma}_f \cdot \mathbf{n}) \cdot \mathbf{t} = \beta \mathbf{u}_f \cdot \mathbf{t} & \text{on } \Gamma_p \\
(\boldsymbol{\sigma}_p \cdot \mathbf{n}) \cdot \mathbf{t} = 0 & \text{on } \Gamma_p \\
\mathbf{u}_f = \mathbf{0} & \text{on } \Gamma_f
\end{array}
\right.
\end{aligned}
\end{equation}
$$

3.3 Stokes-Brinkman coupling with relaxation of interface conditions

We can propose a “regularised” version of the four conditions (23)-(26) by imposing, for a fixed choice of a num-
ber $\varepsilon > 0$, the following interface conditions

$$
\left\{
\begin{array}{ll}
((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{n} = 0 & \text{on } \Gamma_p \\[1em]
((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{t} = \beta(\mathbf{u}_f + \varepsilon^2 \mathbf{u}_p) \cdot \mathbf{t} & \text{on } \Gamma_p \\[1em]
((\boldsymbol{\sigma}_f + \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{n} = \frac{2}{\varepsilon}(\mathbf{u}_f - \mathbf{u}_p) \cdot \mathbf{n} & \text{on } \Gamma_p \\[1em]
((\boldsymbol{\sigma}_f + \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{t} = \beta(\mathbf{u}_f - \varepsilon^2 \mathbf{u}_p) \cdot \mathbf{t} & \text{on } \Gamma_p
\end{array}
\right.
\qquad (28)
$$

The choice of these conditions is due to several reasons.

1. This type of general condition linking strain jumps, velocity jumps as well as stress and velocity averages at the interface can be found in the work of P. Angot [2, 6, 3]. The aim here is to weight the different terms in a clever way.

2. The introduction of the parameter $\varepsilon > 0$ gives more regular solutions, see the theorem 1 and the definitions 1 for $\varepsilon > 0$ compared to the definition 2 corresponding to the case $\varepsilon = 0$. In particular, for all $\varepsilon > 0$, there is no need to impose an interface condition on the test functions in the weak formulation, which simplifies the implementation of a numerical scheme.

3. By performing the sum and difference of equations (29) and (31), we have

$$
(29) \text{ and } (31) \iff \left\{
\begin{array}{ll}
(\boldsymbol{\sigma}_f \cdot \mathbf{n}) \cdot \mathbf{t} = \beta \mathbf{u}_f \cdot \mathbf{t} & \text{on } \Gamma_p \\[1em]
(\boldsymbol{\sigma}_p \cdot \mathbf{n}) \cdot \mathbf{t} = -\beta \varepsilon^2 \mathbf{u}_p \cdot \mathbf{t} & \text{on } \Gamma_p
\end{array}
\right.
$$

so that we can directly observe that we formally find the model (SB) when $\varepsilon$ goes to 0 (a rigorous proof will
be given in Section 4.2).

The complete regularized model writes

$$
\begin{equation}
\begin{split}
\text{(SR)} \quad \left\{
\begin{array}{@{}l@{\quad}l@{}}
-\eta \Delta \mathbf{u}_f + \nabla p_f = \mathbf{F} & \text{in } \Omega_f \\[1ex]
\operatorname{div} \mathbf{u}_f = 0 & \text{in } \Omega_f \\[1ex]
-\gamma \Delta \mathbf{u}_p + \kappa \mathbf{u}_p + \nabla p_p = \mathbf{0} & \text{in } \Omega_p \\[1ex]
\operatorname{div} \mathbf{u}_p = 0 & \text{in } \Omega_p \\[1ex]
((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{n} = 0 & \text{on } \Gamma_p \\[1ex]
((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{t} = \beta (\mathbf{u}_f + \varepsilon^2 \mathbf{u}_p) \cdot \mathbf{t} & \text{on } \Gamma_p \\[1ex]
((\boldsymbol{\sigma}_f + \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{n} = \frac{2}{\varepsilon} (\mathbf{u}_f - \mathbf{u}_p) \cdot \mathbf{n} & \text{on } \Gamma_p \\[1ex]
((\boldsymbol{\sigma}_f + \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{t} = \beta (\mathbf{u}_f - \varepsilon^2 \mathbf{u}_p) \cdot \mathbf{t} & \text{on } \Gamma_p \\[1ex]
\mathbf{u}_f = \mathbf{0} & \text{on } \Gamma_f
\end{array}
\right.
\end{split}
\end{equation}
$$

We now have three models - (SD), (SB), (SR) - coupling Stokes flow in $\Omega_f$ and porous medium flow in $\Omega_p$. The aim of the following section is to propose a mathematical analysis of these models.

4 Mathematical analysis

More precisely, the goal of this section is to show that the three problems (SD), (SB) and (SR) are well posed: existence and uniqueness of solution, in a sense to be defined. We also wish to understand the influence of parameters, typically $\varepsilon$ and $\gamma$, and whether the transition between these models can be made by continuously vanishing these parameters.

The results presented here are a logical continuation of the work done by several teams on relatively similar problems.
In [14], the authors define a weak solution of the coupling of time-dependent incompressible Navier-Stokes equations with Darcy equations where the interface conditions include the Beavers-Joseph-Saffman condition. They prove existence and uniqueness of the weak solution by a constructive approach with weak regularity interfaces. If the fluid flow is parallel to the fluid-porous interface, many authors study the model. As example, in [18], the authors prove the well-posedness of the coupled Stokes-Darcy problem with a new set of effective interface conditions under a suitable relationship between the permeability and the boundary layer constants containing geometrical information about the porous medium and the interface. Recently, in [15], the authors obtain the global existence and uniqueness of the strong solution to an interface problem between a fluid flow, governed by Navier-Stokes equations, and a flow in a porous medium governed by the Darcy equation. Their proof is based on the local well-posedness and some key a priori bounds established. Cao et al in [13] consider the full Beavers-Joseph-Jones conditions and investigate the well-posedness of a coupled Stokes-Darcy model with Beavers-Joseph interface boundary conditions without the convection term. In the steady-state case, the well-posedness is established under the assumption of small coefficient in the Beavers-Joseph interface boundary condition. In the time-dependent case, the well-posedness is established via appropriate time discretization of the problem and a novel scaling of the system under isotropic media assumption.

Following the example of this work, we will first define, for each of the three problems (SD), (SB) and (SR), what we mean by a solution. The key point is, of course, to write down the interface conditions correctly, so that we can hope to show their existence and uniqueness.

## 4.1 Weak formulations

We introduce the functional spaces

$$ 
\begin{aligned}
\mathbf{L}_{\text{div}}^{2}(\Omega_{p}) &= \{\mathbf{v} \in L^{2}(\Omega_{p})^{2} ; \operatorname{div} \mathbf{v} = 0\} \\
\mathbf{H}_{\text{div}}^{1}(\Omega_{p}) &= \{\mathbf{v} \in H^{1}(\Omega_{p})^{2} ; \operatorname{div} \mathbf{v} = 0\} \\
\mathbf{H}_{\text{div}}^{1}(\Omega_{f}) &= \{\mathbf{v} \in H^{1}(\Omega_{f})^{2} ; \operatorname{div} \mathbf{v} = 0 \text{ and } \mathbf{v}|_{\Gamma_{f}} = \mathbf{0}\}
\end{aligned}
$$

equipped with the usual inner product and norm (thanks to the Poincaré's inequality), that is respectively

$$ 
\begin{aligned}
\|\mathbf{v}\|_{0,p}^{2} &= \int_{\Omega_{p}} |\mathbf{v}|^{2}, \\
\|\mathbf{v}\|_{1,p}^{2} &= \int_{\Omega_{p}} |\mathbf{v}|^{2} + \int_{\Omega_{p}} |\mathbf{D}\mathbf{v}|^{2}, \\
\|\mathbf{v}\|_{1,f}^{2} &= \int_{\Omega_{f}} |\mathbf{D}\mathbf{v}|^{2}.
\end{aligned}
$$

We also note

$$ 
\begin{aligned}
\mathbf{W} &= \mathbf{H}_{\text{div}}^{1}(\Omega_{f}) \times \mathbf{H}_{\text{div}}^{1}(\Omega_{p}), \\
\mathbf{W}_{0} &= \{(\mathbf{v}_{f}, \mathbf{v}_{p}) \in \mathbf{W} ; (\mathbf{v}_{f} - \mathbf{v}_{p}) \cdot \mathbf{n} = 0 \text{ on } \Gamma_{p},\} \\
\mathbf{Z}_{0} &= \{(\mathbf{v}_{f}, \mathbf{v}_{p}) \in \mathbf{H}_{\text{div}}^{1}(\Omega_{f}) \times \mathbf{L}_{\text{div}}^{2}(\Omega_{p}) ; (\mathbf{v}_{f} - \mathbf{v}_{p}) \cdot \mathbf{n} = 0 \text{ on } \Gamma_{p}\}
\end{aligned}
$$

that we will equip with the usual product norms.

**Definition 1** Let $\eta > 0, \kappa > 0, \beta > 0, \gamma > 0, \varepsilon > 0$ and $\mathbf{F} \in L^2(\Omega_f)^2$. We say that $(\mathbf{u}_f, \mathbf{u}_p) \in \mathbf{W}$ is a weak solution to the problem (SR) if for all $(\varphi_f, \varphi_p) \in \mathbf{W}$ we have

$$ \mathbf{a}_{\text{SR}}((\mathbf{u}_f, \mathbf{u}_p), (\varphi_f, \varphi_p)) = \mathbf{b}(\varphi_f, \varphi_p) $$

where the bilinear form $\mathbf{a}_{\text{SR}}$ is defined by

$$ 
\begin{aligned}
\mathbf{a}_{\text{SR}}((\mathbf{u}_f, \mathbf{u}_p), (\varphi_f, \varphi_p)) ={}& 2\eta \int_{\Omega_f} \mathbf{D}\mathbf{u}_f : \mathbf{D}\varphi_f + 2\gamma \int_{\Omega_p} \mathbf{D}\mathbf{u}_p : \mathbf{D}\varphi_p + \kappa \int_{\Omega_p} \mathbf{u}_p \cdot \varphi_p \\
& + \frac{1}{\varepsilon} \int_{\Gamma_p} ((\mathbf{u}_f - \mathbf{u}_p) \cdot \mathbf{n})((\varphi_f - \varphi_p) \cdot \mathbf{n}) + \beta \int_{\Gamma_p} (\mathbf{u}_f \cdot \mathbf{t})(\varphi_f \cdot \mathbf{t}) + \beta \varepsilon^2 \int_{\Gamma_p} (\mathbf{u}_p \cdot \mathbf{t})(\varphi_p \cdot \mathbf{t})
\end{aligned}
\quad (41) $$

and where the linear form $\mathbf{b}$ is defined by

$$ \mathbf{b}(\varphi_f, \varphi_p) = \int_{\Omega_f} \mathbf{F} \cdot \varphi_f. \quad (42) $$

**Remark 4** All the terms in $\mathbf{a}_{\text{SR}}$ make sense since any function of $H^1(\Omega_f)$ or $H^1(\Omega_p)$ has a trace on the boundary $\Gamma_p$ which belongs to $L^2(\Gamma_p)$ (and even to $H^{1/2}(\Gamma_p)$ as long as $\Gamma_p$ is Lipschitz).

The definition 1 is justified because all strong solution to (SR) is a weak solution. Indeed, consider $(\mathbf{u}_f, p_f, \mathbf{u}_p, p_p)$ a regular solution of (SR). We perform the scalar product of Equation (32) by $\varphi_f \in$
$\mathbf{H}_{\text{div}}^{1}(\Omega_{f})$ and integrate on $\Omega_{f}$. After integration by parts (remember that $\mathbf{n}$ denotes the unit normal vector at the interface $\Gamma_{p}$, outside the domain $\Omega_{p}$), we obtain

$$2\eta \int_{\Omega_f} \mathbf{D}u_f : \mathbf{D}\varphi_f + \int_{\Gamma_p} (\sigma_f \cdot \mathbf{n}) \cdot \varphi_f = \int_{\Omega_f} \mathbf{F} \cdot \varphi_f. \quad (43)$$

In the same way, we perform the scalar product of Equation (34) by $\varphi_p \in \mathbf{H}_{\text{div}}^{1}(\Omega_p)$ and we integrate on $\Omega_p$:

$$2\gamma \int_{\Omega_p} \mathbf{D}u_p : \mathbf{D}\varphi_p + \kappa \int_{\Omega_p} \mathbf{u}_p \cdot \varphi_p - \int_{\Gamma_p} (\sigma_p \cdot \mathbf{n}) \cdot \varphi_p = 0. \quad (44)$$

The sum (43) and (44) gives the quantity

$$X = (\sigma_f \cdot \mathbf{n}) \cdot \varphi_f - (\sigma_p \cdot \mathbf{n}) \cdot \varphi_p$$

who writes

$$X = \frac{1}{2}((\sigma_f + \sigma_p) \cdot \mathbf{n}) \cdot (\varphi_f - \varphi_p) + \frac{1}{2}((\sigma_f - \sigma_p) \cdot \mathbf{n}) \cdot (\varphi_f + \varphi_p). \quad (45)$$

With the boundary conditions (36)–(39), the equality (45) writes

$$X = \frac{1}{\varepsilon}((\mathbf{u}_f - \mathbf{u}_p) \cdot \mathbf{n})((\varphi_f - \varphi_p) \cdot \mathbf{n}) + \beta(\mathbf{u}_f \cdot \mathbf{t})(\varphi_f \cdot \mathbf{t}) + \beta\varepsilon^2(\mathbf{u}_p \cdot \mathbf{t})(\varphi_p \cdot \mathbf{t}). \quad (46)$$

Then the sum (43)+(44) gives the equality

$$2\eta \int_{\Omega_f} \mathbf{D}u_f : \mathbf{D}\varphi_f + 2\gamma \int_{\Omega_p} \mathbf{D}u_p : \mathbf{D}\varphi_p + \kappa \int_{\Omega_p} \mathbf{u}_p \cdot \varphi_p \\ + \frac{1}{\varepsilon} \int_{\Gamma_p} ((\mathbf{u}_f - \mathbf{u}_p) \cdot \mathbf{n})((\varphi_f - \varphi_p) \cdot \mathbf{n}) + \beta \int_{\Gamma_p} (\mathbf{u}_f \cdot \mathbf{t})(\varphi_f \cdot \mathbf{t}) + \beta\varepsilon^2 \int_{\Gamma_p} (\mathbf{u}_p \cdot \mathbf{t})(\varphi_p \cdot \mathbf{t}) = \int_{\Omega_f} \mathbf{F} \cdot \varphi_f$$

that corresponds to the formulation introduced in the definition 1.

In the same way, we introduce the notion of a weak solution for problem (SB).

**Definition 2** Let $\eta > 0$, $\kappa > 0$, $\beta > 0$, $\gamma > 0$ and $\mathbf{F} \in L^2(\Omega_f)^2$.

We say that $(\mathbf{u}_f, \mathbf{u}_p) \in \mathbf{W}_0$ is a weak solution to the problem (SB) if for all $(\varphi_f, \varphi_p) \in \mathbf{W}_0$ we have

$$\mathbf{a}_{\text{SB}}((\mathbf{u}_f, \mathbf{u}_p), (\varphi_f, \varphi_p)) = \mathbf{b}(\varphi_f, \varphi_p)$$

where the bilinear form $\mathbf{a}_{\text{SR}}$ is defined by

$$\mathbf{a}_{\text{SB}}((\mathbf{u}_f, \mathbf{u}_p), (\varphi_f, \varphi_p)) = 2\eta \int_{\Omega_f} \mathbf{D}u_f : \mathbf{D}\varphi_f + 2\gamma \int_{\Omega_p} \mathbf{D}u_p : \mathbf{D}\varphi_p + \kappa \int_{\Omega_p} \mathbf{u}_p \cdot \varphi_p + \beta \int_{\Gamma_p} (\mathbf{u}_f \cdot \mathbf{t})(\varphi_f \cdot \mathbf{t}).$$

The idea for obtaining this definition is similar to that proposed for the definition 1. We do the scalar product of Equation (19) by $\varphi_p \in \mathbf{H}_{\text{div}}^{1}(\Omega_p)$ and we integrate on $\Omega_p$. After integration by parts, we obtain

$$2\eta \int_{\Omega_f} \mathbf{D}u_f : \mathbf{D}\varphi_f + \int_{\Gamma_p} (\sigma_f \cdot \mathbf{n}) \cdot \varphi_f = \int_{\Omega_f} \mathbf{F} \cdot \varphi_f. \quad (47)$$

Similarly, we calculate the scalar product of Equation (21) by $\varphi_p \in \mathbf{H}_{\text{div}}^{1}(\Omega_p)$ and we integrate on $\Omega_p$:

$$2\gamma \int_{\Omega_p} \mathbf{D}u_p : \mathbf{D}\varphi_p + \kappa \int_{\Omega_p} \mathbf{u}_p \cdot \varphi_p - \int_{\Gamma_p} (\sigma_p \cdot \mathbf{n}) \cdot \varphi_p = 0 \quad (48)$$

However, the expression (46) for the quantity $X$ involved in the terms at the boundary $\Gamma_p$ is no longer valid. We use the conditions (23), (26) and $(\varphi_f - \varphi_p) \cdot \mathbf{n} = 0$ on $\Gamma_p$ (because $(\varphi_f, \varphi_p) \in \mathbf{W}_0$) in order to obtain

$$X = \beta(\mathbf{u}_f \cdot \mathbf{t})(\varphi_f \cdot \mathbf{t}) \quad (49)$$

which explains the formulation for the bilinear form $\mathbf{a}_{\text{SB}}$.

Finally, we also introduce the notion of weak solution to the problem (SD)

**Definition 3** Let $\eta > 0$, $\kappa > 0$, $\beta > 0$ and $\mathbf{F} \in L^2(\Omega_f)^2$.

We say that $(\mathbf{u}_f, \mathbf{u}_p) \in \mathbf{Z}_0$ is a weak solution to the problem (SD) if for all $(\varphi_f, \varphi_p) \in \mathbf{Z}_0$ we have

$$\mathbf{a}_{\text{SD}}((\mathbf{u}_f, \mathbf{u}_p), (\varphi_f, \varphi_p)) = \mathbf{b}(\varphi_f, \varphi_p)$$

where the bilinear form $\mathbf{a}_{\text{SR}}$ is defined by

$$\mathbf{a}_{\text{SD}}((\mathbf{u}_f, \mathbf{u}_p), (\varphi_f, \varphi_p)) = 2\eta \int_{\Omega_f} \mathbf{D}u_f : \mathbf{D}\varphi_f + \kappa \int_{\Omega_p} \mathbf{u}_p \cdot \varphi_p + \beta \int_{\Gamma_p} (\mathbf{u}_f \cdot \mathbf{t})(\varphi_f \cdot \mathbf{t}).$$
## 4.2 Results of existence, uniqueness and convergence

**Theorem 1 (Existence and uniqueness)** Let $\eta > 0$, $\kappa > 0$, $\beta > 0$ and $\mathbf{F} \in L^2(\Omega_f)^2$.

1. For all $\varepsilon > 0$ and $\gamma > 0$ there exists a unique weak solution $(\mathbf{u}_f, \mathbf{u}_p)$, in the sense of Definition 1, to problem (SR).

2. For all $\gamma > 0$ there exists a unique weak solution $(\mathbf{u}_f, \mathbf{u}_p)$, in the sense of Definition 2, to problem (SB).

3. There exists a unique weak solution $(\mathbf{u}_f, \mathbf{u}_p)$, in the sense of Definition 3, to problem (SD).

**Proof:** The proof of these results is performed using the Lax-Milgram's theorem.

The continuity of $\mathbf{b}$ and $\mathbf{a}_{SD}$, $\mathbf{a}_{SB}$ and $\mathbf{a}_{SR}$ are given directly. About the fact that $\mathbf{a}_{SD}$, $\mathbf{a}_{SB}$ and $\mathbf{a}_{SR}$ are coercive we proceed as follows. For $(\mathbf{v}_f, \mathbf{v}_p) \in \mathbf{W}$, we have

$$ \mathbf{a}_{SR}((\mathbf{v}_f, \mathbf{v}_p), (\mathbf{v}_f, \mathbf{v}_p)) \geq \mathbf{a}_{SB}((\mathbf{v}_f, \mathbf{v}_p), (\mathbf{v}_f, \mathbf{v}_p)) \geq 2\eta \|\mathbf{v}_f\|_{1,f}^2 + \min(2\gamma, \kappa) \|\mathbf{v}_p\|_{1,p}^2 $$

The right-hand side of this inequality being the square of a norm on $\mathbf{W}$ as well as on $\mathbf{W}_0$ this implies the coercivity of $\mathbf{a}_{SR}$ and that of $\mathbf{a}_{SB}$. In the same way, for the bilinear form $\mathbf{a}_{SD}$, we note that for $(\mathbf{v}_f, \mathbf{v}_p) \in \mathbf{Z}_0$ we have

$$ \mathbf{a}_{SD}((\mathbf{v}_f, \mathbf{v}_p), (\mathbf{v}_f, \mathbf{v}_p)) \geq 2\eta \|\mathbf{v}_f\|_{1,f}^2 + \kappa \|\mathbf{v}_p\|_{0,p}^2 \quad (50) $$

which ensures the coercivity of $\mathbf{a}_{SD}$.

<figure><img src="image_4.png" /></figure>

**Remark 5** The choice of interface conditions is fundamental in the proof of Theorem 1. For example, for the Stokes-Darcy coupling, if we use the Beaver-Joseph's condition instead of the Beaver-Joseph-Saffman's condition, then the weak formulation gives rise to an boundary term $(\beta \int_{\Gamma} \mathbf{u}_p \cdot \boldsymbol{\varphi}_p)$ which is difficult to control. In this case, however, the paper by Hou and Qin [21] provides a positive answer to existence and uniqueness.

**Theorem 2 (Convergence)** Let $\eta > 0$, $\kappa > 0$, $\beta > 0$ and $\mathbf{F} \in L^2(\Omega_f)^2$.

For $\gamma > 0$ and $\varepsilon > 0$, we denote $(\mathbf{u}_f^{\varepsilon,\gamma}, \mathbf{u}_p^{\varepsilon,\gamma})$ the weak solution to problem (SR), $(\mathbf{u}_f^{0,\gamma}, \mathbf{u}_p^{0,\gamma})$ the weak solution to problem (SB) and $(\mathbf{u}_f^{0,0}, \mathbf{u}_p^{0,0})$ the weak solution to problem (SD).

1. The sequence $(\mathbf{u}_f^{\varepsilon,\gamma}, \mathbf{u}_p^{\varepsilon,\gamma})_{\varepsilon>0}$ weakly converges to $(\mathbf{u}_f^{0,\gamma}, \mathbf{u}_p^{0,\gamma})$ in $\mathbf{W}$ when $\varepsilon$ tends to 0.

2. The sequence $(\mathbf{u}_f^{0,\gamma}, \mathbf{u}_p^{0,\gamma})_{\gamma>0}$ weakly converges to $(\mathbf{u}_f^{0,0}, \mathbf{u}_p^{0,0})$ in $\mathbf{Z}_0$ when $\gamma$ tends to 0.

**Proof of the first step** We choose $(\varphi_f, \varphi_p) = (\mathbf{u}_f^{\varepsilon,\gamma}, \mathbf{u}_p^{\varepsilon,\gamma})$ in the definition 1 of the weak solution to the problem (SR). We obtain

$$ \begin{aligned} 2\eta \int_{\Omega_f} |\mathbf{Du}_f^{\varepsilon,\gamma}|^2 + 2\gamma \int_{\Omega_p} |\mathbf{Du}_p^{\varepsilon,\gamma}|^2 + \kappa \int_{\Omega_p} |\mathbf{u}_p^{\varepsilon,\gamma}|^2 & \\ + \frac{1}{\varepsilon} \int_{\Gamma_p} |(\mathbf{u}_f^{\varepsilon,\gamma} - \mathbf{u}_p^{\varepsilon,\gamma}) \cdot \mathbf{n}|^2 + \beta \int_{\Gamma_p} |\mathbf{u}_f^{\varepsilon,\gamma} \cdot \mathbf{t}|^2 + \beta\varepsilon^2 \int_{\Gamma_p} |\mathbf{u}_p^{\varepsilon,\gamma} \cdot \mathbf{t}|^2 = \int_{\Omega_f} \mathbf{F} \cdot \mathbf{u}_f^{\varepsilon,\gamma}. & \end{aligned} \quad (51) $$

Using Cauchy-Schwarz, Young and Poincaré inequalities, we obtain

$$ \int_{\Omega_f} \mathbf{F} \cdot \mathbf{u}_f^{\varepsilon,\gamma} \leq \eta \int_{\Omega_f} |\mathbf{Du}_f^{\varepsilon,\gamma}|^2 + C $$

where the constant $C$ depends only on $\mathbf{F}$, $\Omega_f$ and $\eta$ but not on $\varepsilon$. We deduce then the following estimate

$$ \begin{aligned} \eta \int_{\Omega_f} |\mathbf{Du}_f^{\varepsilon,\gamma}|^2 + 2\gamma \int_{\Omega_p} |\mathbf{Du}_p^{\varepsilon,\gamma}|^2 + \kappa \int_{\Omega_p} |\mathbf{u}_p^{\varepsilon,\gamma}|^2 & \\ + \frac{1}{\varepsilon} \int_{\Gamma_p} |(\mathbf{u}_f^{\varepsilon,\gamma} - \mathbf{u}_p^{\varepsilon,\gamma}) \cdot \mathbf{n}|^2 + \beta \int_{\Gamma_p} |\mathbf{u}_f^{\varepsilon,\gamma} \cdot \mathbf{t}|^2 + \beta\varepsilon^2 \int_{\Gamma_p} |\mathbf{u}_p^{\varepsilon,\gamma} \cdot \mathbf{t}|^2 \leq C. & \end{aligned} \quad (52) $$

We obtain uniform bounds (with respect to $\varepsilon$) on the quantities $\|\mathbf{Du}_f^{\varepsilon,\gamma}\|_{1,f}$ and $\|\mathbf{Du}_p^{\varepsilon,\gamma}\|_{1,p}$. Using the continuity of the trace maps $\mathbf{H}_{\text{div}}^1(\Omega_f) \hookrightarrow L^2(\Gamma)^2$ we also deduce uniform bounds on the quantities $\|\mathbf{u}_f^{\varepsilon,\gamma}\|_{L^2(\Gamma)}$, and similarly on $\|\mathbf{u}_p^{\varepsilon,\gamma}\|_{L^2(\Gamma)}$.

Such bounds allow us to extract some sub-sequences, always indexed by $\varepsilon$, which weakly converge when $\varepsilon$ tends to 0: there exist $(\bar{\mathbf{u}}_f^{0,\gamma}, \bar{\mathbf{u}}_p^{0,\gamma}) \in \mathbf{W}$ such that

$$ \begin{aligned} \mathbf{Du}_f^{\varepsilon,\gamma} \rightarrow \mathbf{D}\bar{\mathbf{u}}_f^{0,\gamma} & \quad L^2(\Omega_f)\text{-weak,} & \mathbf{u}_f^{\varepsilon,\gamma} \rightarrow \bar{\mathbf{u}}_f^{0,\gamma} & \quad L^2(\Gamma)\text{-weak} \\ \mathbf{Du}_p^{\varepsilon,\gamma} \rightarrow \mathbf{D}\bar{\mathbf{u}}_p^{0,\gamma} & \quad L^2(\Omega_p)\text{-weak,} & \mathbf{u}_p^{\varepsilon,\gamma} \rightarrow \bar{\mathbf{u}}_p^{0,\gamma} & \quad L^2(\Gamma)\text{-weak.} \end{aligned} \quad (53) $$

We note that the compactness of $\mathbf{H}_{\text{div}}^1(\Omega_p)$ in $L^2(\Omega_p)^2$ gives

$$ \mathbf{u}_p^{\varepsilon,\gamma} \rightarrow \bar{\mathbf{u}}_p^{0,\gamma} \quad L^2(\Omega_p)\text{-strong.} \quad (54) $$

Moreover the singular term with $\frac{1}{\varepsilon}$ in the estimate (52) implies $(\bar{\mathbf{u}}_f^{0,\gamma} - \bar{\mathbf{u}}_p^{0,\gamma}) \cdot \mathbf{n} = 0$ then $(\bar{\mathbf{u}}_f^{0,\gamma}, \bar{\mathbf{u}}_p^{0,\gamma}) \in \mathbf{W}_0$.

Finally if we choose $(\varphi_f, \varphi_p) \in \mathbf{W}_0$ in the definition 1 of the weak solution to the problem (SR), the singular term with $\frac{1}{\varepsilon}$ disappears (since $(\varphi_f - \varphi_p) \cdot \mathbf{n} = 0$) and we pass to the limit $\varepsilon$ tends to 0 in the all other terms thanks to the convergence (53) and (54). We deduce that the limit $(\bar{\mathbf{u}}_f^{0,\gamma}, \bar{\mathbf{u}}_p^{0,\gamma})$ is the weak solution to the problem (SB). Thanks to the uniqueness of this solution, we have $(\bar{\mathbf{u}}_f^{0,\gamma}, \bar{\mathbf{u}}_p^{0,\gamma}) = (\mathbf{u}_f^{0,\gamma}, \mathbf{u}_p^{0,\gamma})$ which gives the result.

<figure><img src="image_5.png" /></figure>
**Proof of the second step** We choose $(\varphi_f, \varphi_p) = (\mathbf{u}_f^{0,\gamma}, \mathbf{u}_p^{0,\gamma})$ in the definition 2 of the weak solution to the problem (SB). We obtain

$$2\eta \int_{\Omega_f} |\mathbf{D}\mathbf{u}_f^{0,\gamma}|^2 + 2\gamma \int_{\Omega_p} |\mathbf{D}\mathbf{u}_p^{0,\gamma}|^2 + \kappa \int_{\Omega_p} |\mathbf{u}_p^{0,\gamma}|^2 + \beta \int_{\Gamma_p} |\mathbf{u}_f^{0,\gamma} \cdot \mathbf{t}|^2 = \int_{\Omega_f} \mathbf{F} \cdot \mathbf{u}_f^{0,\gamma}. \quad (55)$$

As previously, the Cauchy-Schwarz, Young and Poincaré's inequalities allow us to control the right-hand side so that we obtain the estimate

$$\eta \int_{\Omega_f} |\mathbf{D}\mathbf{u}_f^{0,\gamma}|^2 + 2\gamma \int_{\Omega_p} |\mathbf{D}\mathbf{u}_p^{0,\gamma}|^2 + \eta \int_{\Omega_p} |\mathbf{u}_p^{0,\gamma}|^2 + \beta \int_{\Gamma_p} |\mathbf{u}_f^{0,\gamma} \cdot \mathbf{t}|^2 \le C \quad (56)$$

where C is a constant independent on $\gamma$.

As previously, such bounds allow us to extract some sub-sequences, always indexed by $\gamma$, which converge when $\gamma$ tends to 0: there exists $(\bar{\mathbf{u}}_f^{0,0}, \bar{\mathbf{u}}_p^{0,0}) \in \mathbf{W}_0$ such that

$$\begin{array}{lll}
\mathbf{D}\mathbf{u}_f^{0,\gamma} \to \mathbf{D}\bar{\mathbf{u}}_f^{0,0} & L^2(\Omega_f)\text{-weak,} & \mathbf{u}_f^{0,\gamma} \to \bar{\mathbf{u}}_f^{0,0} \quad L^2(\Gamma)\text{-weak} \\
\gamma\mathbf{D}\mathbf{u}_p^{0,\gamma} \to \mathbf{0} & L^2(\Omega_p)\text{-weak,} & \mathbf{u}_p^{0,\gamma} \to \bar{\mathbf{u}}_p^{0,0} \quad L^2(\Omega_p)\text{-weak}
\end{array}$$

We note that the weak limit $\gamma \mathbf{D}\bar{\mathbf{u}}_p^{0,\gamma} \to \mathbf{0}$ comes from the uniform bound on $\sqrt{\gamma} \mathbf{D}\mathbf{u}_p^{0,\gamma}$.

We pass to the limit $\gamma$ tends to 0 in the weak formulation $a_{SB}((\mathbf{u}_f^{0,\gamma}, \mathbf{u}_p^{0,\gamma}), (\varphi_f, \varphi_p)) = b(\varphi_f, \varphi_p)$. We deduce that $a_{SD}((\bar{\mathbf{u}}_f^{0,0}, \bar{\mathbf{u}}_p^{0,0}), (\varphi_f, \varphi_p)) = b(\varphi_f, \varphi_p)$ and then $(\bar{\mathbf{u}}_f^{0,0}, \bar{\mathbf{u}}_p^{0,0}) = (\mathbf{u}_f^{0,0}, \mathbf{u}_p^{0,0})$ is the unique weak solution to the problem (SD). $\square$

# 5 Taking inertia effects into account

In this section, we show that it is possible to take into account the nonlinear term $\mathbf{u}_f \cdot \nabla \mathbf{u}_f$ in the air flow model. More precisely, we are interested in the Navier-Stokes-Darcy model containing this convection term

$$\text{(NSD)} \quad \left\{
\begin{array}{ll}
\rho \mathbf{u}_f \cdot \nabla \mathbf{u}_f - \eta \Delta \mathbf{u}_f + \nabla p_f = \mathbf{F} & \text{in } \Omega_f \\
\text{div } \mathbf{u}_f = 0 & \text{in } \Omega_f \\
\kappa \mathbf{u}_p + \nabla p_p = \mathbf{0} & \text{in } \Omega_p \\
\text{div } \mathbf{u}_p = 0 & \text{in } \Omega_p \\
(\mathbf{u}_f - \mathbf{u}_p) \cdot \mathbf{n} = 0 & \text{on } \Gamma_p \\
((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{n} = 0 & \text{on } \Gamma_p \\
((\boldsymbol{\sigma}_f - \boldsymbol{\sigma}_p) \cdot \mathbf{n}) \cdot \mathbf{t} = \beta \mathbf{u}_f \cdot \mathbf{t} & \text{on } \Gamma_p \\
\mathbf{u}_f = \mathbf{0} & \text{on } \Gamma_f
\end{array}
\right.$$

Of course, the same analysis can be extended to the Stokes-Brinkman (SB) model or the regularized (SR) model, without further difficulty.

**Definition 4** Let $\rho > 0$, $\eta > 0$, $\kappa > 0$, $\beta > 0$ and $\mathbf{F} \in L^2(\Omega_f)^2$.

We say that $(\mathbf{u}_f, \mathbf{u}_p) \in \mathbf{Z}_0$ is a weak solution to the (NSD) problem if, for all $(\varphi_f, \varphi_p) \in \mathbf{Z}_0$ we have

$$\rho \int_{\Omega_f} (\mathbf{u}_f \cdot \nabla \mathbf{u}_f) \cdot \varphi_f + a_{SD}((\mathbf{u}_f, \mathbf{u}_p), (\varphi_f, \varphi_p)) = b(\varphi_f, \varphi_p)$$

where we recall that the bilinear form $a_{SR}$ is defined by

$$a_{SD}((\mathbf{u}_f, \mathbf{u}_p), (\varphi_f, \varphi_p)) = 2\eta \int_{\Omega_f} \mathbf{D}\mathbf{u}_f : \mathbf{D}\varphi_f + \kappa \int_{\Omega_p} \mathbf{u}_p \cdot \varphi_p + \beta \int_{\Gamma_p} (\mathbf{u}_f \cdot \mathbf{t})(\varphi_f \cdot \mathbf{t}),$$

and that the linear for $b$ is defined by $b(\varphi_f, \varphi_p) = \int_{\Omega_f} \mathbf{F} \cdot \varphi_f$.

The following theorem ensures the existence and the uniqueness of a solution to this problem when the density $\rho$ is not too large. Note that this type of condition seems natural, since it is also required in the case of stationary Navier-Stokes equations. More precisely, the condition we are going to impose involves the injection constant $C_s > 0$ from $H^1(\Omega_f) \subset L^4(\Omega_f)$ defined by

$$\exists C_s > 0; \forall u \in H^1(\Omega_f) \quad \|u\|_{L^4(\Omega_f)} \le C_s \|u\|_{H^1(\Omega_f)},$$

whose value essentially depends on the domain $\Omega_f$.

**Theorem 3 (Existence and uniqueness)** Let $\rho > 0$, $\eta > 0$, $\kappa > 0$, $\beta > 0$ and $\mathbf{F} \in L^2(\Omega_f)^2$.

If $\rho C_s^2 \|\mathbf{F}\|_{L^2(\Omega_f)} < 4\eta^2$ then there exists a unique weak solution $(\mathbf{u}_f, \mathbf{u}_p)$, in the sense of Definition 4, to (NSD) problem.
**Proof:** For $R > 0$, we define $B_R = \{(\mathbf{u}_f, \mathbf{u}_p) \in \mathbb{Z}_0 ; \|\mathbf{u}_f\|_{L^4(\Omega_f)} < R\}$. The value of the real number $R$ will be specified later depending on the requirements. The proof of Theorem 3 will be split into three parts. First we consider a linearised version of (NSD) for which there is a unique solution. We then obtain bounds on this solution that allow us to show, via a fixed point theorem, the announced result.

**Step 1 - linear problem associated to the (NSD) nonlinear problem**

Let $(\bar{\mathbf{u}}_f, \bar{\mathbf{u}}_p) \in B_R$ and consider the problem: find $(\mathbf{u}_f, \mathbf{u}_p) \in \mathbb{Z}_0$ such that for all $(\varphi_f, \varphi_p) \in \mathbb{Z}_0$ we have

$$ \rho \int_{\Omega_f} (\bar{\mathbf{u}}_f \cdot \nabla \mathbf{u}_f) \cdot \varphi_f + \mathbf{a}_{SD}((\mathbf{u}_f, \mathbf{u}_p), (\varphi_f, \varphi_p)) = \mathbf{b}(\varphi_f, \varphi_p). \quad (57) $$

The existence of a unique solution for such a problem comes from to the Lax-Milgram's theorem. Indeed, we will prove that the bilinear application $\mathbf{a}_{NSD}: \mathbb{Z}_0 \times \mathbb{Z}_0 \to \mathbb{R}$ defined by

$$ \mathbf{a}_{NSD}((\mathbf{u}_f, \mathbf{u}_p), (\varphi_f, \varphi_p)) = \rho \int_{\Omega_f} (\bar{\mathbf{u}}_f \cdot \nabla \mathbf{u}_f) \cdot \varphi_f + \mathbf{a}_{SD}((\mathbf{u}_f, \mathbf{u}_p), (\varphi_f, \varphi_p)) $$

is continuous and coercive.

✓ The application $\mathbf{a}_{NSD}$ is continuous because we already know that $\mathbf{a}_{SD}$ is continuous and because

$$ \begin{aligned} \left| \rho \int_{\Omega_f} (\bar{\mathbf{u}}_f \cdot \nabla \mathbf{u}_f) \cdot \varphi_f \right| & \le \rho \|\bar{\mathbf{u}}_f\|_{L^4(\Omega_f)} \|\nabla \mathbf{u}_f\|_{L^2(\Omega_f)} \|\varphi_f\|_{L^4(\Omega_f)} \\ & \le \rho C_s R \|\mathbf{u}_f\|_{1,f} \|\varphi_f\|_{1,f}. \end{aligned} $$

✓ Using the same kind of estimate, we prove that the application $\mathbf{a}_{NSD}$ is coercive since

$$ \left| \rho \int_{\Omega_f} (\bar{\mathbf{u}}_f \cdot \nabla u_f) \cdot \mathbf{u}_f \right| \le \rho C_s \|\bar{\mathbf{u}}_f\|_{L^4} \|\mathbf{u}_f\|_{1,f}^2 $$

which, thanks to Inequality (50), implies

$$ \begin{aligned} \mathbf{a}_{NSD}((\mathbf{u}_f, \mathbf{u}_p), (\mathbf{u}_f, \mathbf{u}_p)) & \ge (2\eta - \rho C_s \|\bar{\mathbf{u}}_f\|_{L^4}) \|\mathbf{u}_f\|_{1,f}^2 + \kappa \|\mathbf{u}_p\|_{0,p}^2 \\ & \ge (2\eta - \rho C_s R) \|\mathbf{u}_f\|_{1,f}^2 + \kappa \|\mathbf{u}_p\|_{0,p}^2. \end{aligned} $$

Finally $\mathbf{a}_{NSD}$ is coercive as soon as $2\eta - \rho C_s R > 0$. From now, we assume that $R$ is small enough in order to satisfy such inequality.

**Step 2 - bound on the solution of the linear problem**

Using $(\varphi_f, \varphi_p) = (\mathbf{u}_f, \mathbf{u}_p)$ in the weak formulation (57), we obtain the following energy estimate on the solution

$$ \begin{aligned} 2\eta \|\mathbf{u}_f\|_{1,f}^2 + \kappa \|\mathbf{u}_p\|_{0,p}^2 & \le \|\mathbf{F}\|_{L^2(\Omega_f)} \|\mathbf{u}_f\|_{1,f} + \rho \|\bar{\mathbf{u}}_f\|_{L^4} \|\nabla \mathbf{u}_f\|_{L^2} \|\mathbf{u}_f\|_{L^4} \\ & \le \|\mathbf{F}\|_{L^2(\Omega_f)} \|\mathbf{u}_f\|_{1,f} + \rho RC_s \|\mathbf{u}_f\|_{1,f}^2. \end{aligned} $$

We deduce in particular that

$$ \|\mathbf{u}_f\|_{1,\Omega} \le \frac{1}{2\eta - \rho RC_s} \|\mathbf{F}\|_{L^2(\Omega_f)}. \quad (58) $$

**Step 3 - nonlinear problem solution as a fixed point**

Now, using the first step, we define the mapping

$$ \Phi(\bar{\mathbf{u}}_f, \bar{\mathbf{u}}_p) = (\mathbf{u}_f, \mathbf{u}_p). $$

Let us check that, for $R$ small enough, $\Phi$ is a contraction on $B_R$. Consider $(\bar{\mathbf{u}}_f, \bar{\mathbf{u}}_p) \in B_R$, $(\bar{\mathbf{v}}_f, \bar{\mathbf{v}}_p) \in B_R$, denote $(\mathbf{u}_f, \mathbf{u}_p) = \Phi(\bar{\mathbf{u}}_f, \bar{\mathbf{u}}_p)$, $(\mathbf{v}_f, \mathbf{v}_p) = \Phi(\bar{\mathbf{v}}_f, \bar{\mathbf{v}}_p)$ and consider the differences

$$ (\mathbf{w}_f, \mathbf{w}_p) = (\mathbf{u}_f - \mathbf{v}_f, \mathbf{u}_p - \mathbf{v}_p) $$

and

$$ (\bar{\mathbf{w}}_f, \bar{\mathbf{w}}_p) = (\bar{\mathbf{u}}_f - \bar{\mathbf{v}}_f, \bar{\mathbf{u}}_p - \bar{\mathbf{v}}_p). $$

The goal is then to prove that there exists $K < 1$, independent of $\mathbf{w}_f, \mathbf{w}_p, \bar{\mathbf{w}}_f$ and $\bar{\mathbf{w}}_p$, such that

$$ \| (\mathbf{w}_f, \mathbf{w}_p) \|_{\mathbb{Z}_0} \le K \|(\bar{\mathbf{w}}_f, \bar{\mathbf{w}}_p)\|_{\mathbb{Z}_0} $$

where the norm on $\mathbb{Z}_0$ will be explicitly precised hereafter.

By making the difference between the weak formulation satisfied by $(\mathbf{u}_f, \mathbf{u}_p)$ and by $(\mathbf{v}_f, \mathbf{v}_p)$, we deduce that for all $(\varphi_f, \varphi_p) \in \mathbb{Z}_0$, we have

$$ \rho \int_{\Omega} (\bar{\mathbf{w}}_f \cdot \nabla \mathbf{u}_f) \cdot \varphi_f + \rho \int_{\Omega} (\bar{\mathbf{v}}_f \cdot \nabla \mathbf{w}_f) \cdot \varphi_f + \mathbf{a}_{SD}((\mathbf{w}_f, \mathbf{w}_p), (\varphi_f, \varphi_p)) = 0. $$
We choose $\varphi_f = \mathbf{w}_f$ and $\varphi_p = \mathbf{w}_p$. Due to the inequality (50) and the Hölder inequality, we have

$$2\eta \|\mathbf{w}_f\|_{1,f}^2 + \kappa \|\mathbf{w}_p\|_{0,p}^2 \le \rho \|\bar{\mathbf{w}}_f\|_{L^4} \|\nabla \mathbf{u}_f\|_{L^2} \|\mathbf{w}_f\|_{L^4} + \rho \|\bar{\mathbf{v}}_f\|_{L^4} \|\nabla \mathbf{w}_f\|_{L^2} \|\mathbf{w}_f\|_{L^4}.$$

Using the continuous injection $H^1(\Omega_f) \subset L^4(\Omega_f)$ again, and the fact that $\|\bar{\mathbf{v}}_f\|_{L^4} < R$ we deduce

$$(2\eta - \rho C_s R) \|\mathbf{w}_f\|_{1,f}^2 + \kappa \|\mathbf{w}_p\|_{0,p}^2 \le \rho C_s^2 \|\bar{\mathbf{w}}_f\|_{1,f} \|\mathbf{u}_f\|_{1,f} \|\mathbf{w}_f\|_{1,f}.$$

Using the bound (58), we obtain

$$(2\eta - \rho C_s R) \|\mathbf{w}_f\|_{1,f}^2 + \kappa \|\mathbf{w}_p\|_{0,p}^2 \le \frac{\rho C_s^2}{2\eta - \rho RC_s} \|\mathbf{F}\|_{L^2(\Omega_f)} \|\bar{\mathbf{w}}_f\|_{1,f} \|\mathbf{w}_f\|_{1,f} \quad (59)$$

Particularly, we deduce

$$\|\mathbf{w}_f\|_{1,f} \le \frac{\rho C_s^2}{(2\eta - \rho RC_s)^2} \|\mathbf{F}\|_{L^2(\Omega_f)} \|\bar{\mathbf{w}}_f\|_{1,f},$$

so that Estimation (59) becomes

$$(2\eta - \rho C_s R) \|\mathbf{w}_f\|_{1,f}^2 + \kappa \|\mathbf{w}_p\|_{0,p}^2 \le \frac{\rho^2 C_s^4}{(2\eta - \rho RC_s)^3} \|\mathbf{F}\|_{L^2(\Omega_f)}^2 \|\bar{\mathbf{w}}_f\|_{1,f}^2. \quad (60)$$

We can rewrite this estimate as follows

$$\|(\mathbf{w}_f, \mathbf{w}_p)\|_{\mathbf{Z}_0}^2 \le \frac{\rho^2 C_s^4}{(2\eta - \rho RC_s)^4} \|\mathbf{F}\|_{L^2(\Omega_f)}^2 \|(\bar{\mathbf{w}}_f, \bar{\mathbf{w}}_p)\|_{\mathbf{Z}_0}^2 \quad (61)$$

where we use the following norm on $\mathbf{Z}_0$

$$\|(\mathbf{w}_f, \mathbf{w}_p)\|_{\mathbf{Z}_0}^2 = (2\eta - \rho C_s R) \|\mathbf{w}_f\|_{1,f}^2 + \kappa \|\mathbf{w}_p\|_{0,p}^2.$$

Consequently, the mapping $\Phi$ is a contraction on $B_R$ if

$$\rho C_s^2 \|\mathbf{F}\|_{L^2(\Omega_f)} < (2\eta - \rho RC_s)^2$$

which is satisfied for $R$ small enough as soon as $\rho C_s^2 \|\mathbf{F}\|_{L^2(\Omega_f)} < 4\eta^2$ exactly corresponding to the announced condition. Therefore, $\Phi$ has exactly one fixed point $(\mathbf{u}_f, \mathbf{u}_p) \in \mathbf{Z}_0$ which is the unique weak solution of (NSD). $\square$

# 6 Numerical simulations and comments

In this section, we give some simulations. We have numerical experimented in an urban environment, in order to assess the impacts. More precisely, we consider a rectangular domain containing two sub-domains modelling buildings (in black on the figure 5). The air flow velocity at the building boundaries is assumed to be zero (homogeneous Dirichlet condition), while the velocity at the outer boundaries is assumed to be equal to (1, 0) and corresponds to an imposed horizontal velocity (see Figure 2, left).

<figure><img src="image_6.png" /><figcaption>Figure 2: Left: Geometrical configuration for simulations: two buildings in grey and wind imposed on the boundary. Right: Velocity norm for the treeless case and wind intensity scale.</figcaption></figure>

As a reference, we simulated an air flow with this geometry and Figure 2 (right) shows the result in terms of wind intensity. This intensity scale will be used for all subsequent simulations.

**Remark 6** All the simulations were carried out using the Freefem++ software, in which the weak form of the equations was implemented. Post-processing is performed using ParaView software.
<table><caption>Table 1: Values of parameter *k* as a function of *d* and *φ*.</caption><thead><tr><th>d</th><th>0.8</th><th>0.9</th><th>0.99</th></tr></thead><tbody><tr><th>0.01</th><td>8.5 · 10<sup>-6</sup></td><td>4.9 · 10<sup>-5</sup></td><td>6.5 · 10<sup>-3</sup></td></tr><tr><th>0.05</th><td>2.1 · 10<sup>-4</sup></td><td>1.2 · 10<sup>-3</sup></td><td>1.6 · 10<sup>-1</sup></td></tr><tr><th>0.1</th><td>8.5 · 10<sup>-4</sup></td><td>4.9 · 10<sup>-3</sup></td><td>6.5 · 10<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 3: Configurations with different value for the permeability parameter *k*. From top to bottom, left to right: *k* = 1.6 · 10⁻¹, *k* = 5 · 10⁻³, *k* = 3 · 10⁻³, *k* = 2 · 10⁻³, *k* = 1.2 · 10⁻³.</figcaption></figure>

## 6.1 Influence of the permeability wood

As previously presented, the Kozeny-Carman's law (7) expresses the permeability *k* with respect to the length *d* and the porosity *φ*, quantities which both have a precise physical meaning. Indeed, an analogy can be made with the permeability of a homogeneous granular medium, where the length *d* corresponds to the grain diameter and the quantity *φ* to the porosity. The following table gives values that seem realistic for these quantities, as well as the value for the corresponding coefficient *k*.

We therefore tested different values of this parameter *k* in a given geometric configuration (and *ε* = 0.01). Figure 3 shows the configuration and the results according to the different values of *k*. We note in particular that the smaller the parameter *k*, the more the forest acts as a brake on the propagation of velocity within it.

## 6.2 Influence of the ε parameter

In this section, we tested different values of the parameter *ε* for a given geometric configuration and permeability. In future work, we will enrich these models by taking into account the thermal and evapotranspiration effects of wooded areas. The results, given in figure 4, show that this parameter plays much the same role as permeability.

## 6.3 Influence of the geometrical configuration of wooded areas

Finally, we wanted to compare different scenarios for the establishment of wooded areas. The three scenarios are shown in figure 5 while the results (for the parameters *k* = 1.2 · 10⁻³ and *ε* = 0.01) are given in figure 6.
<figure><img src="image_8.png" /><figcaption>Figure 4: Configurations with different value for the parameter $\epsilon$. From top to bottom, left to right: $\epsilon = 0.006$, $\epsilon = 0.01$, $\epsilon = 0.02$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 5: Configuration without plantation (top/left : two buildings in black, wind imposed on the boundary) and three planting scenarios</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 6: Results for the four scenarios</figcaption></figure>

## References

[1] F.M. Allan and M.H. Hamdan. Fluid mechanics of the interface region between two porous layers. *Applied Mathematics and computation*, 128(1):37-43, 2002.
[2] P. Angot. On the well-posed coupling between free fluid and porous viscous flows. *Applied mathematics letters*, 24(6):803–810, 2011.

[3] P. Angot. Well-posed Stokes/Brinkman and Stokes/Darcy coupling revisited with new jump interface conditions. *ESAIM: Mathematical Modelling and Numerical Analysis*, 52(5):1875–1911, 2018.

[4] P. Angot, C.H. Bruneau, and P. Fabrie. A penalization method to take into account obstacles in incompressible viscous flows. *Numerische Mathematik*, 81(4):497–520, 1999.

[5] P. Angot, G. Carbou, and V. Péron. Asymptotic study for Stokes-Brinkman model with jump embedded transmission conditions. *Asymptotic Analysis*, 96(3-4):223–249, 2016.

[6] P. Angot, B. Goyeau, and J.A. Ochoa-Tapia. Asymptotic modeling of transport phenomena at the interface between a fluid and a porous layer: Jump conditions. *Physical Review E*, 95(6):063302, 2017.

[7] P. Angot, B. Goyeau, and J.A. Ochoa-Tapia. A nonlinear asymptotic model for the inertial flow at a fluid-porous interface. *Advances in Water Resources*, 149:103798, 2021.

[8] A. Bastide, P.H. Cocquet, and D. Ramalingam. Penalization model for Navier-Stokes-Darcy equations with application to porosity-oriented topology optimization. *Mathematical Models and Methods in Applied Sciences*, 28(08):1481–1512, 2018.

[9] G. Beavers and D.D. Joseph. Boundary conditions at a naturally permeable wall. *Journal of fluid mechanics*, 30(1):197–207, 1967.

[10] M.E. Bogovskii. Solution of the first boundary value problem for the equation of continuity of an incompressible medium. In *Doklady Akademii Nauk*, volume 5-248, pages 1037–1040. Russian Academy of Sciences, 1979.

[11] F. Boyer and P. Fabrie. *Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models*, volume 183. Springer Science & Business Media, 2012.

[12] C.H. Bruneau and P. Fabrie. New efficient boundary conditions for incompressible Navier-Stokes equations: a well posedness result. *M2AN*, 30(7):815–840, 1996.

[13] Y. Cao, M. Gunzburger, F. Hua, and X. Wang. Coupled Stokes-Darcy model with Beavers-Joseph interface boundary condition. *Commun. Math. Sci.*, 8-1:1–25, 2010.

[14] A. Cesmelioglu, V. Girault, and B. Riviere. Time-dependent coupling of Navier-Stokes and Darcy flows. *ESAIM: Mathematical Modelling and Numerical Analysis*, 47(2):539–554, 2013.

[15] M. Cui, W. Dong, and Z. Guo. Global well-posedness of coupled Navier-Stokes and Darcy equations. *Journal of Differential Equations*, 388:82–111, 2024.

[16] M. Discacciati, A. Quarteroni, et al. Navier-Stokes/Darcy coupling: modeling, analysis, and numerical approximation. *Rev. Mat. Complut*, 22(2):315–426, 2009.

[17] L. Durlofsky and J.F. Brady. Analysis of the Brinkman equation as a model for flow in porous media. *Physics of Fluids*, 30(11):3329, 1987.

[18] E. Eggenweiler, M. Discacciati, and I. Rybak. Analysis of the Stokes-Darcy problem with generalised interface conditions. *ESAIM: Mathematical Modelling and Numerical Analysis*, 56(2):727–742, 2022.

[19] H.I. Ene et al. Equations et phénomènes de surface pour l'écoulement dans un modèle de milieu poreux. *J. Mec.*, 14-1:73–108, 1975.

[20] J. Fuchsberger, P. Aigner, S. Niederer, G. Plank, H. Schima, G. Haase, and E. Karabelas. On the incorporation of obstacles in a fluid flow problem using a Navier-Stokes-Brinkman penalization approach. *Journal of Computational Science*, 57:101506, 2022.

[21] Y. Hou and Y. Qin. On the solution of coupled stokes/darcy model with Beavers-Joseph interface condition. *Computers & Mathematics with Applications*, 77(1):50–65, 2019.

[22] W.J. Layton, F. Schieweck, and I. Yotov. Coupling fluid flow with porous media flow. *SIAM Journal on Numerical Analysis*, 40(6):2195–2218, 2002.

[23] A. Mikelic and W. Jäger. On the interface boundary condition of Beavers, Joseph, and Saffman. *SIAM Journal on Applied Mathematics*, 60(4):1111–1127, 2000.

[24] D.A. Nield. The Beavers-Joseph boundary condition and related matters: a historical and critical note. *Transport in porous media*, 78:537–540, 2009.

[25] D.A. Nield, A. Bejan, et al. *Convection in porous media*, volume 3. Springer, 2006.

[26] M. Parvazinia, V. Nassehi, R.J. Wakeman, and M.H.R. Ghoreishy. Finite element modelling of flow through a porous medium between two parallel plates using the Brinkman equation. *Transport in porous media*, 63:71–90, 2006.

[27] P.G. Saffman. On the boundary condition at the surface of a porous medium. *Studies in applied mathematics*, 50(2):93–101, 1971.

[28] P. Strohbeck, E. Eggenweiler, and I. Rybak. A modification of the Beavers-Joseph condition for arbitrary flows to the fluid-porous interface. *Transport in Porous Media*, 147(3):605–628, 2023.