1 Introduction

Coupled porous-medium and free-flow systems appear in a variety of environmental and technical applications such as interactions between surface and subsurface water, geothermal problems and industrial filtration processes (Beaude et al. 2019; Cimolin and Discacciati 2013; Reuter et al. 2019). From the macroscale perspective, fluid flow in the free-flow domain and through the porous medium is described by two different sets of equations, which need to be coupled at the fluid–porous interface (Angot et al. 2017; Discacciati et al. 2002; Jäger and Mikelić 2009; Lācis and Bagheri 2017; Beavers and Joseph 1967). Fluid flows through such coupled systems are highly interface driven, therefore, the correct choice of coupling conditions and effective model parameters is crucial for physically consistent modeling and accurate numerical simulation of applications.

Depending on the flow problem, different mathematical models and different sets of coupling conditions at the fluid–porous interface are applied, e.g., (Dawson 2008; Angot et al. 2017; Jäger and Mikelić 2009; Lācis and Bagheri 2017; Carraro et al. 2015; Sochala et al. 2009; Eggenweiler and Rybak 2020; Discacciati and Quarteroni 2009; Girault and Rivière 2009). Most often, the Navier–Stokes equations or, in case of low Reynolds numbers, the Stokes equations are used to describe fluid flow in the free-flow region, while in the porous-medium domain, Darcy’s law or one of its generalizations is considered. Besides the Stokes equations, there are several other simplifications of the Navier–Stokes equations, which can be used to model surface flows, such as the shallow water equations, kinematic or diffusive wave equations (Magiera et al. 2016; Sochala et al. 2009; Dawson 2008). For higher Reynolds numbers, the Forchheimer equation is used in the porous medium and coupled to the free-flow model (Angot et al. 2017; Cimolin and Discacciati 2013). To describe two-fluid-phase flows in the subsurface, the Richards equation or two-phase Darcy’s law is applied (Mosthaf et al. 2011; Beaude et al. 2019; Rybak et al. 2015).

Based on the choice of the underlying flow models in the free-flow region and in the porous medium, different sets of coupling conditions at the fluid–porous interface are available in the literature. Interface conditions for the Navier–Stokes/Darcy–Forchheimer model are developed in Cimolin and Discacciati (2013); Amara et al. (2009); Alazmi and Vafai (2001); Angot et al. (2021). Coupling conditions for the Stokes equations and multiphase Darcy’s law together with the transport of chemical species and energy are proposed in Mosthaf et al. (2011). Coupling strategies for the Stokes/Darcy–Brinkman equations are derived in Angot et al. (2017). The most widely studied problem in the literature, both from the modeling and numerical point of view, is the Stokes/Darcy problem describing single-fluid-phase flows in coupled free-flow and porous-medium systems. For this model, different coupling strategies exist, e.g., (Beavers and Joseph 1967; Angot et al. 2017; Lācis and Bagheri 2017; Jäger and Mikelić 2000; Eggenweiler and Rybak 2021; Ochoa-Tapia and Whitaker 1995). However, the classical set of coupling conditions, which comprises the conservation of mass across the interface, the balance of normal forces and the Beavers–Joseph condition on the tangential velocity (Beavers and Joseph 1967), is routinely used. Often, the Saffman simplification of the Beavers–Joseph condition neglecting the porous-medium velocity at the interface is considered (Saffman 1971). An extension of the Beavers–Joseph coupling condition with the symmetrized viscous stress tensor was postulated by Jones (1973).

Originally, the Beavers–Joseph coupling condition was proposed for flows parallel to the fluid–porous interface. Nevertheless, this condition and also its variants obtained by Saffman (1971) and Jones (1973) are often used for non-parallel flows to the porous layer, e.g., (Hanspal et al. 2006; Discacciati and Gerardo-Giorda 2018). In Eggenweiler and Rybak (2020), it is shown that the Beavers–Joseph condition is not applicable for arbitrary flow directions to the interface that is the case for, e.g., filtration problems. Besides the unsuitability of the Beavers–Joseph coupling condition (10) for non-parallel flows to the interface, the slip coefficient \(\alpha _{\textrm{BJ}}\) is uncertain and needs to be determined for every flow problem. Although the question concerning the correct value of the Beavers–Joseph parameter is known for decades, to the best of the authors knowledge, no systematic study on determination of this parameter has been published. In case of cylindrical grains and parallel flows to the porous bed, the values of the Beavers–Joseph slip coefficient were obtained in Mierzwiczak et al. (2019). In the literature, \(\alpha _{\textrm{BJ}}=1\) is typically used even if this is not the correct choice (Yang et al. 2019; Eggenweiler and Rybak 2020; Rybak et al. 2021; Lācis et al. 2020). However, the correct value of the Beavers–Joseph parameter is essential for accurate numerical simulations of applications, where the Beavers–Joseph condition is applicable, e.g., for microfluidic experiments (Terzis et al. 2019).

Alternative coupling concepts for the Stokes/Darcy problem exist in the literature. However, some of them contain unknown model parameters, which still need to be determined before the interface conditions can be used in numerical simulations of applications, e.g., (Angot et al. 2017, 2021; Ochoa-Tapia and Whitaker 1995). Other alternative coupling concepts, where the effective coefficients can be computed based on the pore geometry, are either derived only for unidirectional flows parallel or perpendicular to the porous layer (Jäger and Mikelić 2000, 2009; Carraro et al. 2015) or they are not validated for arbitrary flow directions (Lācis and Bagheri 2017; Lācis et al. 2020; Zampogna and Bottaro 2016). An exception is the set of generalized interface conditions proposed in Eggenweiler and Rybak (2021), which is valid for arbitrary flows in coupled free-flow and porous-medium systems. This set of coupling conditions comprises the conservation of mass across the interface (8) or (11), the classical balance of normal forces (9) for isotropic porous media and its extension (12) in case of anisotropic media, and a generalization (13) of the Beavers–Joseph condition (10) on the tangential velocity. Condition (13) contains an interfacial porous-medium velocity instead of the Darcy velocity, which appears in the Beavers–Joseph condition (10), and the boundary layer coefficient \({\varvec{\mathrm N}}^{\textrm{bl}}\) instead of the Beavers–Joseph slip parameter \(\alpha _{\textrm{BJ}}\). The goal of this paper is to reformulate the generalized interface condition (13) in the same analytical form as the traditionally used Beavers–Joseph condition (10) and to provide the means to compute effective coefficients appearing in the modified condition. In this way, other researchers can use their already developed software for coupled flow problems based on the classical interface conditions and only adjust the effective model parameters.

In this paper, we study two different flow scenarios: a pressure-driven flow parallel to the porous layer and a general filtration problem with arbitrary flow directions to the interface. We formulate the generalized interface condition (13) on the tangential velocity component in the form of the Beavers–Joseph condition (10) and determine the effective model parameters for various porous-medium geometrical configurations (isotropic, orthotropic, anisotropic). To compute the effective coefficients for periodic media, the theory of homogenization and boundary layers (Jäger and Mikelić 1996; Hornung 1997) is applied. As an alternative approach, we propose an efficient two-level numerical algorithm, which is valid also for non-periodic porous structures. We demonstrate that the Beavers–Joseph slip coefficient \(\alpha _{\textrm{BJ}}\) can be found only for unidirectional flows to the interface and compute the optimal value of this parameter for all considered pore geometries. We show that \(\alpha _{\textrm{BJ}}\) cannot be fitted for arbitrary flow directions, and in this case, modification (14) has to be used. We perform model validation by comparison of pore-scale resolved to macroscale numerical simulation results.

The paper is organized as follows. In Section 2, we describe the pore-scale resolved and the macroscale flow models including two sets of coupling conditions on the fluid–porous interface: the classical and the modified ones. In Section 3, we provide the approach to compute permeability and boundary layer constants appearing in the interface conditions and propose an efficient two-level numerical algorithm to determine the optimal value of the Beavers–Joseph slip coefficient. The detailed numerical study of different coupled flow problems is conducted in Section 4 and the conclusions are given in Section 5.

2 Mathematical Models

In this section, we describe the flow system of interest and present the microscale (pore-scale resolved) and the macroscale flow models. The microscale model is used for the determination of the optimal Beavers–Joseph slip coefficient for unidirectional flows, for the computation of the effective coefficients in the modified conditions, and for the validation of interface conditions. The macroscale model consists of the Stokes/Darcy equations with the classical set of coupling conditions and the modified ones, respectively.

2.1 Assumptions on Geometry and Flow

We consider steady-state fluid flow in the coupled domain \(\Omega = \Omega _{\textrm{ff}}\cup \Omega _{\textrm{pm}}\) consisting of the free-flow region \(\Omega _{\textrm{ff}}\subset \mathbb {R}^2\) and the porous medium \(\Omega _{\text {pm}} \subset \mathbb {R}^2\). The porous layer contains periodically distributed solid inclusions that allows to compute the permeability and the effective boundary layer coefficients in the modified interface conditions using the theory of homogenization. We assume the separation of scales, i.e., \(\ell \ll L\), where \(\ell\) denotes the characteristic pore size and L is the length of the domain (Fig. 1).

Fig. 1
figure 1

Flow system description at the microscale (left) and macroscale (right)

Full size image

We consider incompressible single-fluid-phase flows at low Reynolds numbers (\(Re \ll 1\)). The porous medium is supposed to be fully saturated with the same fluid that is present in the free-flow region. From the microscale perspective, the Stokes equations describe the fluid flow in the whole flow domain \(\Omega _{\textrm{mi}}=\Omega _{\textrm{ff}}\cup \Omega _{\textrm{ps}}\), which comprises the free-flow region \(\Omega _{\textrm{ff}}\) and the pore space \(\Omega _{\textrm{ps}}\) of the porous medium (Fig. 1, left). At the macroscale, the flow system is divided into two different continuum flow regions \(\Omega _{\textrm{ff}}\) and \(\Omega _{\textrm{pm}}\) separated by the sharp fluid–porous interface \(\Gamma\) (Fig. 1, right), which is void of thermodynamic properties. In this case, two different models are applied in the two flow domains and the models need to be coupled on the interface \(\Gamma\).

2.2 Microscale Model

Under the assumptions given in Section 2.1, fluid flow in the whole flow domain \(\Omega _{\textrm{mi}}\) is governed by the Stokes equations

$$\begin{aligned} \nabla \varvec{\mathsf {\cdot }}{\varvec{\mathrm v}}_{\textrm{mi}}= 0, \qquad - \nabla \varvec{\mathsf {\cdot }}\varvec{\textsf{T}}({\varvec{\mathrm v}}_{\textrm{mi}},p_{\textrm{mi}}) - \rho {\varvec{\mathrm g}} = {\varvec{\mathrm 0}} \qquad \text {in} \;\; \Omega _{\textrm{mi}}, \end{aligned}$$
(1)

completed with the no-slip condition on the fluid–solid boundary

$$\begin{aligned} {\varvec{\mathrm v}}_{\textrm{mi}}= {\varvec{\mathrm 0}} \quad \text {on} \;\; \partial \Omega _{\textrm{mi}}\setminus \partial \Omega , \end{aligned}$$
(2)

and appropriate conditions on the external boundary \(\partial \Omega = \Gamma _D \cup \Gamma _N\):

$$\begin{aligned} {\varvec{\mathrm v}}_{\textrm{mi}}= \overline{{\varvec{\mathrm v}}} \quad \text {on} \;\; \Gamma _{D}, \qquad \varvec{\textsf{T}}({\varvec{\mathrm v}}_{\textrm{mi}}, \, p_{\textrm{mi}}) {\varvec{\mathrm n}} = \overline{{\varvec{\mathrm h}}} \quad \text {on} \;\; \Gamma _{N}. \end{aligned}$$
(3)

Here, \({\varvec{\mathrm v}}_{\textrm{mi}}\) and \(p_{\textrm{mi}}\) denote the fluid velocity and pressure, \(\rho\) is the fluid density, \({\varvec{\mathrm g}}\) is the gravitational acceleration, \(\varvec{\textsf{T}}({\varvec{\mathrm v}},p) = \mu \nabla {\varvec{\mathrm v}} - p {\varvec{\textsf{I}}}\) is the stress tensor, \({\varvec{\textsf{I}}}\) is the identity tensor, \(\mu\) is the dynamic viscosity, \(\overline{{\varvec{\mathrm v}}}\) and \(\overline{{\varvec{\mathrm h}}}\) are given functions, and \({\varvec{\mathrm n}}\) is the unit outward normal vector on \(\partial \Omega\).

2.3 Macroscale Model

At the macroscale, we consider two different flow models in the two domains. Fluid flow in the free-flow region \(\Omega _{\textrm{ff}}\) is described by the Stokes equations

$$\begin{aligned} \nabla \varvec{\mathsf {\cdot }}{\varvec{\mathrm v}}_{\textrm{ff}}= 0, \qquad - \nabla \varvec{\mathsf {\cdot }}\varvec{\textsf{T}}({\varvec{\mathrm v}}_{\textrm{ff}},p_{\textrm{ff}}) - \rho {\varvec{\mathrm g}} = {\varvec{\mathrm 0}} \qquad \text {in} \;\; \Omega _{\textrm{ff}}, \end{aligned}$$
(4)

where \({\varvec{\mathrm v}}_{\textrm{ff}}\) is the fluid velocity and \(p_{\textrm{ff}}\) is the pressure. On the external boundary of the free-flow region \(\partial \Omega _{\textrm{ff}}{\setminus }\Gamma = \Gamma _{D, {\textrm{ff}}} \cup \Gamma _{N, {\textrm{ff}}}\), we impose the following boundary conditions

$$\begin{aligned} \begin{aligned} {\varvec{\mathrm v}}_{\textrm{ff}}= \overline{{\varvec{\mathrm v}}}\quad \text {on} \;\; \Gamma _{D, {\textrm{ff}}}, \qquad \varvec{\textsf{T}}({\varvec{\mathrm v}}_{\textrm{ff}}, \, p_{\textrm{ff}}) {\varvec{\mathrm n}}_{\textrm{ff}}= \overline{{\varvec{\mathrm h}}} \quad \text {on} \;\; \Gamma _{N, {\textrm{ff}}}, \end{aligned} \end{aligned}$$
(5)

where \({\varvec{\mathrm n}}_{\textrm{ff}}\) is the unit outward normal vector on \(\partial \Omega _{\textrm{ff}}\).

Fluid flow through the porous medium \(\Omega _{\textrm{pm}}\) is described by the Darcy flow equations

$$\begin{aligned} \nabla \varvec{\mathsf {\cdot }}{\varvec{\mathrm v}}_{\textrm{pm}}= 0, \qquad {\varvec{\mathrm v}}_{\textrm{pm}}= -\frac{\varvec{\textsf{K}}}{\mu } \left( \nabla p_{\textrm{pm}}- \rho {\varvec{\mathrm g}} \right) \qquad \text {in} \;\; \Omega _{\textrm{pm}}, \end{aligned}$$
(6)

where \({\varvec{\mathrm v}}_{\textrm{pm}}\) is the Darcy velocity, \(p_{\textrm{pm}}\) is the fluid pressure and \(\varvec{\textsf{K}}\) is the intrinsic permeability tensor, which is symmetric, positive definite and bounded. The following boundary conditions are prescribed on the external boundary of the porous-medium domain \(\partial \Omega _{\textrm{pm}}{\setminus } \Gamma = \Gamma _{D, {\textrm{pm}}} \cup \Gamma _{N, {\textrm{pm}}}\):

$$\begin{aligned} \begin{aligned} p_{\textrm{pm}}= {\overline{p}} \quad \text {on} \;\; \Gamma _{D, {\textrm{pm}}}, \qquad {\varvec{\mathrm v}}_{\textrm{pm}}\varvec{\mathsf {\cdot }}{\varvec{\mathrm n}}_{\textrm{pm}}= {\overline{v}} \quad \text {on} \;\; \Gamma _{N, {\textrm{pm}}}, \end{aligned} \end{aligned}$$
(7)

where \(\Gamma _{D, {\textrm{pm}}} \cap \Gamma _{N, {\textrm{pm}}}=\emptyset\), \({\varvec{\mathrm n}}_{\textrm{pm}}\) is the unit outward normal vector on \(\partial \Omega _{\textrm{pm}}\), and \({\overline{p}}\), \({\overline{v}}\) are given functions. Within this work, we neglect the gravitational effects setting \({\varvec{\mathrm g}} = {\varvec{\mathrm 0}}\) in equations (1), (4) and (6).

In addition to the boundary conditions on the external boundary of the coupled domain \(\partial \Omega\), appropriate interface conditions have to be set on the fluid–porous interface \(\Gamma\). In this paper, we consider two different sets of coupling conditions: i) classical conditions which are valid for parallel flows to the interface and ii) generalized conditions which are applicable for arbitrary flow directions.

The classical set of coupling conditions which is most often used in the literature consists of the conservation of mass across the interface (8), the balance of normal forces (9) and the Beavers–Joseph condition (10) for the tangential velocity component

$$\begin{aligned} {\varvec{\mathrm v}}_{\textrm{ff}}\varvec{\mathsf {\cdot }}{\varvec{\mathrm n}}&= {\varvec{\mathrm v}}_{\textrm{pm}}\varvec{\mathsf {\cdot }}{\varvec{\mathrm n}} \qquad{} & {} \text {on} \;\Gamma \, , \end{aligned}$$
(8)
$$\begin{aligned} -{\varvec{\mathrm n}} \varvec{\mathsf {\cdot }}\varvec{\textsf{T}}\left( {\varvec{\mathrm v}}_{\textrm{ff}}, p_{\textrm{ff}} \right) {\varvec{\mathrm n}}&= p_{\textrm{pm}}\qquad{} & {} \text {on} \; \Gamma \, , \end{aligned}$$
(9)
$$\begin{aligned} \left( {\varvec{\mathrm v}}_{\textrm{ff}}-{\varvec{\mathrm v}}_{\textrm{pm}} \right) \varvec{\mathsf {\cdot }}{\varvec{\tau }} - \frac{\sqrt{K}}{\alpha _{\textrm{BJ}}} \nabla {\varvec{\mathrm v}}_{\textrm{ff}}{\varvec{\mathrm n}} \varvec{\mathsf {\cdot }}{\varvec{\tau }}&= 0 \qquad{} & {} \text {on} \; \Gamma \, . \end{aligned}$$
(10)

Here, \({\varvec{\mathrm n}}= -{\varvec{\mathrm n}}_{\textrm{ff}}={\varvec{\mathrm n}}_{\textrm{pm}}\) is the unit vector normal to the fluid–porous interface \(\Gamma\) pointing outward from the porous-medium domain \(\Omega _{\textrm{pm}}\), \({\varvec{\tau }}\) is the unit vector tangential to the interface (Fig. 1, right), \(\alpha _{\textrm{BJ}}>0\) is the Beavers–Joseph slip coefficient and \(K= {\varvec{\tau }} \varvec{\mathsf {\cdot }}\varvec{\textsf{K}} {\varvec{\tau }} \in {\mathbb {R}}^+\).

The generalized interface conditions are recently developed in Eggenweiler and Rybak (2021) for arbitrary flow directions to the fluid–porous interface using homogenization and boundary layer theory. They read

$$\begin{aligned} {\varvec{\mathrm v}}_{\textrm{ff}}\varvec{\mathsf {\cdot }}{\varvec{\mathrm n}}&= {\varvec{\mathrm v}}_{\textrm{pm}}\varvec{\mathsf {\cdot }}{\varvec{\mathrm n}} \qquad{} & {} \text {on} \;\Gamma \, , \end{aligned}$$
(11)
$$\begin{aligned} p_{\textrm{pm}}&= -{\varvec{\mathrm n}} \varvec{\mathsf {\cdot }}\varvec{\textsf{T}}\left( {\varvec{\mathrm v}}_{\textrm{ff}}, p_{\textrm{ff}} \right) {\varvec{\mathrm n}} - \mu N_s^{{\textrm{bl}}} \nabla {\varvec{\mathrm v}}_{\textrm{ff}}{\varvec{\mathrm n}} \varvec{\mathsf {\cdot }}{\varvec{\tau }} \qquad{} & {} \text {on} \;\Gamma \, , \end{aligned}$$
(12)
$$\begin{aligned} {\varvec{\mathrm v}}_{\textrm{ff}}\varvec{\mathsf {\cdot }}{\varvec{\tau }}&= \ell ({\varvec{\mathrm N}}^{{\textrm{bl}}} \varvec{\mathsf {\cdot }}{\varvec{\tau }}) \nabla {\varvec{\mathrm v}}_{\textrm{ff}}{\varvec{\mathrm n}} \varvec{\mathsf {\cdot }}{\varvec{\tau }} - \mu ^{-1}\ell ^2 \sum _{j=1}^2 \left( {\varvec{\mathrm M}}^{j,{\textrm{bl}}} \varvec{\mathsf {\cdot }}{\varvec{\tau }} \right) \frac{\partial p_{\textrm{pm}}}{\partial x_j} \qquad{} & {} \text {on} \;\Gamma \, , \end{aligned}$$
(13)

where \(N_s^{{\textrm{bl}}} \in {\mathbb {R}}, \, {\varvec{\mathrm N}}^{{\textrm{bl}}} = (N_1^{{\textrm{bl}}}, N_2^{{\textrm{bl}}})^\top \in {\mathbb {R}}^2\) and \({\varvec{\mathrm M}}^{j,{\textrm{bl}}} = (M_1^{j,{\textrm{bl}}}, M_2^{j,{\textrm{bl}}})^\top \in {\mathbb {R}}^{2}\) are the boundary layer coefficients with \(N_1^{{\textrm{bl}}}>0\) and \(M_1^{1,{\textrm{bl}}}>0\).

The interface condition (11) is the conservation of mass across the fluid–porous interface, coupling condition (12) is an extension to the balance of normal forces (9), and condition (13) is a generalization of the Beavers–Joseph condition (10) on the tangential velocity component. In this paper, we reformulate condition (13) such that it has the same analytical form as the original Beavers–Joseph interface condition (10):

$$\begin{aligned} \left( {\varvec{\mathrm v}}_{\textrm{ff}}-{\varvec{\mathrm v}}_{\textrm{pm}}^\text {int} \right) \varvec{\mathsf {\cdot }}{\varvec{\tau }} - \ell ({\varvec{\mathrm N}}^{{\textrm{bl}}} \varvec{\mathsf {\cdot }}{\varvec{\tau }}) \nabla {\varvec{\mathrm v}}_{\textrm{ff}}{\varvec{\mathrm n}} \varvec{\mathsf {\cdot }}{\varvec{\tau }} = 0 \qquad \text {on} \;\Gamma \, . \end{aligned}$$
(14)

Here, we define the interfacial velocity as

$$\begin{aligned} {\varvec{\mathrm v}}_{\textrm{pm}}^\text {int} = - \frac{\ell ^2 \varvec{\textsf{M}}^{{\textrm{bl}}}}{\mu } \nabla p_{\textrm{pm}}\, , \end{aligned}$$
(15)

where \(\varvec{\textsf{M}}^{{\textrm{bl}}} \in {\mathbb {R}}^{2 \times 2}\) can be interpreted as the interfacial permeability tensor

$$ \begin{aligned} \varvec{\textsf{M}}^{{\textrm{bl}}} = \begin{pmatrix} M_1^{1,{\textrm{bl}}} \, &{} M_1^{2,{\textrm{bl}}} \\ {} \\M_2^{1,{\textrm{bl}}} \, &{} M_2^{2,{\textrm{bl}}} \end{pmatrix} \, . \end{aligned}$$

Remark 1

In Sudhakar et al. (2021), alternative coupling conditions for the Stokes/Darcy problem are developed. There are similarities between these conditions and the generalized interface conditions (11), (12) and (13) or (14), respectively. The coupling condition on the tangential velocity from Sudhakar et al. (2021) includes an interfacial porous-medium velocity instead of the Darcy velocity as the generalized condition (14). The conservation of mass across the interface (11) is extended to allow transport along \(\Gamma\). The third coupling condition from Sudhakar et al. (2021) is the generalized condition (12) with two additional terms. One term accounts for the normal force induced at the interface due to wall-normal velocity and the other one is a higher-order term that arises during the derivation procedure. The interface conditions developed by Sudhakar et al. (2021) are supposed to account for arbitrary flow directions to the interface, but so far, they have only been validated for the lid-driven cavity over porous bed, where the fluid flow is almost parallel to the interface.

Remark 2

The boundary layer constants \(N_s^{{\textrm{bl}}}, \, {\varvec{\mathrm N}}^{{\textrm{bl}}}\) and \({\varvec{\mathrm M}}^{j,{\textrm{bl}}}\) appearing in the generalized conditions (12) and (13) are computed based on the pore geometry and exact location of the fluid–porous interface as described in Sect. 3. For the particular choice of the interface position, we obtain the corresponding boundary layer constants, which will change if another interface position is considered. Note that, the sharp interface location can be chosen with some freedom above the solid obstacles (Eggenweiler and Rybak 2021).

Remark 3

For orthotropic porous media with \(\varvec{\textsf{K}} = {\text {diag}}\{k_{11}, k_{22}\}\), we have \(N_s^{\textrm{bl}}=0\). In this case, condition (12) reduces to the classical balance of normal forces across the fluid–porous interface given in (9).

Remark 4

The coupling condition (14) contains the interfacial velocity \({\varvec{\mathrm v}}^\text {int}_{\textrm{pm}}\), which is typically different from the Darcy velocity \({\varvec{\mathrm v}}_{\textrm{pm}}\) appearing in the Beavers–Joseph condition (10). For some porous-medium configurations, it is possible to obtain \({\varvec{\mathrm v}}_{\textrm{pm}}^\text {int} \varvec{\mathsf {\cdot }}{\varvec{\tau }} = {\varvec{\mathrm v}}_{\textrm{pm}}\varvec{\mathsf {\cdot }}{\varvec{\tau }}\) by the appropriate choice of the fluid–porous interface location while computing the boundary layer constants \({\varvec{\mathrm M}}^{j,{\textrm{bl}}}\). However, this is not the case in general.

Remark 5

For some orthotropic porous media, which satisfy the condition \({\varvec{\mathrm v}}_{\textrm{pm}}^\text {int} \varvec{\mathsf {\cdot }}{\varvec{\tau }} = {\varvec{\mathrm v}}_{\textrm{pm}}\varvec{\mathsf {\cdot }}{\varvec{\tau }}\) (see Remark 4), we can compute the Beavers–Joseph slip coefficient by setting \(\alpha _{\textrm{BJ}}= \sqrt{K}/(\ell N_1^{\textrm{bl}})\). In this case, the generalized interface condition (14) represents the original Beavers–Joseph condition (10).

3 Computation of Effective Coefficients

In order to obtain computable macroscale models, all effective model parameters need to be determined. For the Stokes/Darcy model (4)–(7) with the classical interface conditions (8)–(10), these are the permeability tensor \(\varvec{\textsf{K}}\) and the Beavers–Joseph slip coefficient \(\alpha _{\textrm{BJ}}\). The permeability tensor is computed using homogenization theory (Hornung 1997) based on the pore geometry. The Beavers–Joseph slip coefficient \(\alpha _{\textrm{BJ}}\) appearing in condition (10) is an uncertain parameter, which needs to be found experimentally or fitted. This parameter depends on the exact position of the fluid–porous interface and the pore geometry in the vicinity of the interface. In this paper, we consider the interface located tangent to the top of the first row of solid inclusions as suggested in Beavers and Joseph (1967); Rybak et al. (2021); Lācis and Bagheri (2017).

Effective coefficients appearing in the Stokes/Darcy model (4)–(7) with the generalized coupling conditions (11)–(13) are the permeability \(\varvec{\textsf{K}}\) and the boundary layer constants \(N_s^{{\textrm{bl}}}, \, {\varvec{\mathrm N}}^{{\textrm{bl}}}\) and \({\varvec{\mathrm M}}^{j,{\textrm{bl}}}\). All these coefficients are computed using homogenization and boundary layer theory (Eggenweiler and Rybak 2021) as described in Sect. 3.1.

In case of fluid flows parallel to the porous bed, the Beavers–Joseph coupling condition (10) is valid. For such flow regimes, we propose an efficient two-level algorithm to determine the optimal value of the Beavers–Joseph coefficient in Sect. 3.2. Taking the optimal value \(\alpha _{\textrm{BJ}}^\text {opt}\) instead of the most commonly used \(\alpha _{\textrm{BJ}}=1\) significantly reduces the error between the pore-scale resolved and the macroscale solutions (see Sect. 4).

Fig. 2
figure 2

Unit cell Y (left) and stripe \(Z^{\textrm{bl}}= Z^+\cup Z^-\) for the solution of the cell problems (17) and the boundary layer problems (18)–(20) and (18)–(19), (22), respectively

Full size image

3.1 Homogenization and Boundary Layer Theory

In order to compute the permeability tensor \(\varvec{\textsf{K}}= (k_{ij})_{i,j=1,2}\), we apply the theory of homogenization (Hornung 1997). This leads to

$$\begin{aligned} k_{ij} =\ell ^2 {\tilde{k}}_{ij}, \qquad {\tilde{k}}_{ij}= \int _{Y_\text {f}} w_{i}^j \ \text {d} {\varvec{\mathrm y}}, \end{aligned}$$
(16)

where \({\varvec{\mathrm w}}^{j} = (w_1^j,w_2^j)\) is the solution to the cell problem

$$\begin{aligned} \begin{aligned} \nabla \varvec{\mathsf {\cdot }}\;&{\varvec{\mathrm w}}^{j} =0 \,, \quad - \Delta {\varvec{\mathrm w}}^{j} + \nabla \pi ^j = {\varvec{\mathrm e}}_j \quad \text {in }Y_\text {f}, \quad \int _{Y_\text {f}} \pi ^j \ \text {d} {\varvec{\mathrm y}} =0 \,, \\&{\varvec{\mathrm w}}^{j} = {\varvec{\mathrm 0}} \quad \text {on }\partial Y_\text {f} \setminus \partial Y \,, \quad \{ {\varvec{\mathrm w}}^{j}, \pi ^j \} \text { is 1-periodic} \,, \quad \end{aligned} \end{aligned}$$
(17)

for \(j=1,2\). Here, \(Y_\text {f}\) denotes the fluid part of the periodicity cell \(Y=(0,1)\times (0,1) = Y_\text {f} \cup Y_\text {s}\) (Fig. 2, left).

The effective coefficients \(N_s^{{\textrm{bl}}}\), \({\varvec{\mathrm N}}^{{\textrm{bl}}}\) and \({\varvec{\mathrm M}}^{j,{\textrm{bl}}}\) appearing in conditions (12) and (13) are obtained from the boundary layer problems constructed in Eggenweiler and Rybak (2021). These problems are defined on the stripe \(Z^{\textrm{bl}}= Z^+\cup Z^-\), where the free-flow part \(Z^+=(0,1)\times (s,4)\) and the pore space \(Z^-= \displaystyle \cup _{k=1}^4 (Y_\text {f}-(0,k))\) are separated by the sharp interface \(S=(0,1) \times \{s \}\) with \(s < 0\). In this paper, we locate the interface S directly on top of the solid inclusions (Fig. 2, right).

We denote the solid boundary within the boundary layer stripe by \(\partial Z^-_\text {s} = \partial Z^- {\setminus } (\partial Z^{\textrm{bl}}\cup S)\), the upper and lower boundary by \(\Gamma _m\) and \(\Gamma _k\), respectively. To obtain coefficients \(N_s^{{\textrm{bl}}}\) and \({\varvec{\mathrm N}}^{{\textrm{bl}}}\), we solve the boundary layer problem

$$\begin{aligned}&\nabla \varvec{\mathsf {\cdot }}{\varvec{\mathrm t}}^{{\textrm{bl}}} = 0 \, , \ \ -\Delta {\varvec{\mathrm t}}^{{\textrm{bl}}} + \nabla s^{{\textrm{bl}}} = {\varvec{\mathrm 0}} \quad \text {in } Z^+\cup Z^-\, , \quad \int _{\Gamma _k} s^{{\textrm{bl}}} \, \text {d}y_1 =0 \, , \quad \end{aligned}$$
(18)
$$\begin{aligned}&\frac{\partial t_1^{{\textrm{bl}}}}{\partial y_2} = t_2^{{\textrm{bl}}} = 0 \quad \text {on } \Gamma _m \, , \quad {\varvec{\mathrm t}}^{{\textrm{bl}}} = {\varvec{\mathrm 0}} \quad \text {on } \partial Z^-_\text {s}\, , \quad \{{\varvec{\mathrm t}}^{{\textrm{bl}}}, s^{{\textrm{bl}}}\} \text { is }y_1- \text {periodic} \, , \end{aligned}$$
(19)
$$\begin{aligned}&\llbracket ( \nabla {\varvec{\mathrm t}}^{{\textrm{bl}}} - s^{{\textrm{bl}}} \varvec{\textsf{I}} ) {\varvec{\mathrm e}}_2\rrbracket _{S} = -{\varvec{\mathrm e}}_1 \, , \ \ \llbracket {\varvec{\mathrm t}}^{{\textrm{bl}}} \rrbracket _{S} = {\varvec{\mathrm 0}} \quad \text {on } S \, , \quad \ {\varvec{\mathrm t}}^{{\textrm{bl}}} = {\varvec{\mathrm 0}} \quad \text {on } \Gamma _k \, , \end{aligned}$$
(20)

and integrate the solutions afterward

$$\begin{aligned} {\varvec{\mathrm N}}^{{\textrm{bl}}}=\left( \int _S t_1^{{\textrm{bl}}} \ \text {d} y_1, 0 \right) \, , \qquad N_s^{\textrm{bl}}= \int _S s^{{\textrm{bl}}} \ \text {d} y_1\, . \end{aligned}$$
(21)

The boundary layer constants \({\varvec{\mathrm M}}^{j,{\textrm{bl}}}\) are obtained by solving two flow problems (18)–(19) with the following jump conditions across the interface S:

$$\begin{aligned} \llbracket ( \nabla {\varvec{\mathrm t}}^{{\textrm{bl}}} - s^{{\textrm{bl}}} \varvec{\textsf{I}} ) {\varvec{\mathrm e}}_2\rrbracket _{S} = \left( \nabla {\varvec{\mathrm w}}^{j} - \pi ^j\varvec{\textsf{I}} \right) {\varvec{\mathrm e}}_2 \, , \quad \llbracket {\varvec{\mathrm t}}^{{\textrm{bl}}} \rrbracket _{S} = {\varvec{\mathrm w}}^j -{{\tilde{k}}}_{2j}{\varvec{\mathrm e}}_2 \ \, \text {on } S \, , \end{aligned}$$
(22)

and integrating their solutions

$$\begin{aligned} {\varvec{\mathrm M}}^{j,{\textrm{bl}}}=\left( \int _S t_1^{{\textrm{bl}}} \ \text {d} y_1, 0 \right) \, , \quad j=1,2 \, . \end{aligned}$$
(23)

3.2 Two-level Numerical Algorithm

In this section, we propose an efficient two-level numerical algorithm to determine the optimal value of the Beavers–Joseph slip coefficient \(\alpha _{\textrm{BJ}}\) for coupled flow problems, where the Beavers–Joseph condition (10) is valid. The algorithm is based on the minimization of the difference between the pore-scale resolved and the macroscale numerical simulation results. This difference is quantified by the errors \(\epsilon _{f, \, c, \, \alpha }\) defined along the fixed cross section \(x_1=c\) as follows

$$\begin{aligned} \begin{aligned} \epsilon _{f, \, c, \, \alpha } = \frac{\Vert f_{{\textrm{ma}}, \, \alpha }(c, \cdot )-f_{\textrm{mi}}(c, \cdot )\Vert }{\Vert f_{\textrm{mi}}(c, \cdot )\Vert }, \qquad \Vert f(c, \cdot )\Vert ^2=\int _h^H f(c,x_2)^2 \, \textrm{d}x_2 , \end{aligned} \end{aligned}$$
(24)

where \(f_{\textrm{mi}}\) is the solution of the pore-scale resolved problem (1)–(3) and \(f_{{\textrm{ma}}, \, \alpha }\) is the solution of the macroscale problem (4)–(10) with \(\alpha _{\textrm{BJ}}= \alpha\). In case the generalized conditions (11)–(13) are applied, we consider \(f_{{\textrm{ma}}}\) for the macroscale solution and denote the error by \(\epsilon _{f, \, c}\).

The Clough–Tocher interpolation is applied to handle the two different meshes used for numerical simulations of pore-scale resolved and macroscale models. Due to the presence of solid obstacles at the microscale, an adaptive triangular mesh is applied in the whole fluid domain \(\Omega _{\textrm{mi}}\) for the pore-scale resolved simulations. The macroscale solutions are computed on the staggered Cartesian mesh in the coupled domain \(\Omega\). To compare the microscale and macroscale solutions, we interpolate the pore-scale simulation results using the Clough–Tocher interpolation (Alfeld 1984; Nielson 1983), which provides \(C^1\)-functions and operates on triangulated data.

Instead of solving the Stokes/Darcy problem (4)–(10) for a huge range of parameters \(\alpha _{\textrm{BJ}}\) to determine the one yielding the minimal error (24), we propose the following efficient two-level procedure. In the first step (Level 1) of the two-level algorithm, we compute a coarse approximation of the optimal value of the Beavers–Joseph parameter. This value is determined in the second refinement step (Level 2).

Level 1: We perform macroscale numerical simulations for the Stokes/Darcy problem (4)–(10) considering the following range of the Beavers–Joseph slip coefficient \({\mathcal {A}}_{\textrm{BJ}}= \{0.01, 1, 2,\ldots , M\}\) with \(M \in {\mathbb {N}}\) and compute the corresponding relative errors (24). We consider the cubic spline interpolation of the \((M+1)\) error values and apply simulated annealing (Alg. 1) to find a coarse approximation \(\alpha _{\textrm{BJ}}^*\) of the optimal Beavers–Joseph parameter. This value is not accurate enough due to interpolation, therefore, the refinement step (Level 2) is applied. We set the lower boundary \(\alpha _{\textrm{BJ}}= 0.01\) due to positivity of the Beavers–Joseph parameter \(\alpha _{\textrm{BJ}}>0\).

figure a

In Algorithm 1, we choose the following monotonically decreasing null sequence \(T_n = a{\text {exp}}(-n/b)\) with \(a=0.005\) and \(b=100\). We set \(c_{\text {grow}} = 1.1\), \(c_{\text {shrink}}=0.99\), \(\varepsilon _\text {tol} = 10^{-12}\) and \(N_\text {max}=10^5\).

Level 2: If the value of the Beavers–Joseph parameter determined in Level 1 is \(\alpha _{\textrm{BJ}}^* \in [0.6,M]\), we round \(\alpha _{\textrm{BJ}}^*\) to the first digit after comma, set \(m=10\) and consider the following candidates for the optimal slip coefficient \({\text {round}}(\alpha _{\textrm{BJ}}^*,1) + 0.5 - j/m\). Otherwise, we consider the candidates 0.01 and \(\lceil \alpha _{\textrm{BJ}}^*\rceil - j/m\) for \(j= 0, \ldots , m-1\). The optimal value \(\alpha _{\textrm{BJ}}^{\text {opt}}\) of the Beavers–Joseph parameter is the candidate which yields the smallest relative error (24). Since the error is changing only slightly when considering more candidates (\(m>10\)) for the numerical simulation results presented in Section 4, it is sufficient to determine \(\alpha _{\textrm{BJ}}^{\text {opt}}\) up to the first digit after comma.

figure b

Remark 6

The upper boundary M for the Beavers–Joseph coefficient \(\alpha _{\textrm{BJ}}\) should be chosen depending on the application of interest. The number of candidates m in Level 2 can be increased in order to determine \(\alpha _{\textrm{BJ}}^\text {opt}\) more accurately, if necessary. Alternative interpolation techniques can be considered in Level 1.

4 Numerical Results

In this section, we study the applicability of the classical interface conditions (8)–(10) and the generalized interface conditions (11)–(13), where condition (13) is formulated in the same form (14) as the original Beavers–Joseph condition (10), to describe fluid flows in coupled free-flow and porous-medium systems. We consider two flow problems: i) pressure-driven flow, where fluid flow is parallel to the porous layer, and ii) general filtration, where flow is arbitrary to the fluid–porous interface. We show that one can find an optimal value of the Beavers–Joseph slip coefficient \(\alpha _{\textrm{BJ}}\) appearing in condition (10) for unidirectional flows parallel to the interface. We determine this value using the efficient numerical algorithm proposed in Sect. 3.2. Alternatively, one can also apply the generalized coupling conditions (11), (12) and (14), where all model parameters are computed based on the pore geometry using homogenization and boundary layer theory (see Sect. 3.1). We demonstrate that for arbitrary flows to the porous layer, it is not possible to find an optimal value of the Beavers–Joseph coefficient. In this case, the coupling condition  (14) should be used instead of the classical Beavers–Joseph condition (10).

We study various porous-medium geometrical configurations considering different shapes of solid inclusions (circular, rectangular, elliptical) and different arrangements of solid grains (in-line, staggered) leading to isotropic, orthotropic and anisotropic media. We consider different porous-medium geometries, which lead to the same permeability and similar porosity (Tab. 1, geometries \(G_1\) and \(G_3\)); porous media containing the same type of solid inclusions of different sizes which result in different interface roughness, permeability and porosity (Tab. 1, geometries \(G_1\), \(G_2\) and \(G_4\)); and pore geometries having the same porosity and the same interface roughness but different permeability tensors (Tab. 1, geometries \(G_1\), \(G_4\) and \(G_5\), \(G_6\)). All considered porous-medium geometries are periodic such that homogenization theory is applied to compute the permeability tensor and the boundary layer constants as described in Sect. 3.1. We set the dynamic viscosity of the fluid \(\mu =1\) for all flow problems.

Table 1 Permeability values \(k_{ij}\), porosity \(\phi\) and boundary layer constants for different pore geometries
Full size table

The microscale problem (1)–(3), the unit cell Stokes problems (17) and the boundary layer problems (18)–(20) and (18), (19), (22) are solved using the software package FreeFEM++ with the Taylor–Hood finite elements (Hecht 2012). The macroscale Stokes/Darcy problem (4)–(10) completed by the classical interface conditions (4)–(7) or the generalized conditions (11), (12), (14) is discretized using the finite volume method on staggered Cartesian grids conforming at the fluid–porous interface with the grid size \(h_1=h_2=1/800\). The optimal value of the Beavers–Joseph slip coefficient \(\alpha _{\textrm{BJ}}^\text {opt}\) is computed using the proposed two-level numerical algorithm (Algorithm 2).

4.1 Pressure-driven Flow

In this section, we study a flow scenario where the flow is parallel to the fluid–porous interface (pressure-driven flow). We consider the free-flow domain \(\Omega _{\textrm{ff}}=[0,1]\times [0,0.5]\), the porous medium \(\Omega _{\textrm{pm}}=[0,1]\times [-0.5,0]\) and the sharp fluid–porous interface \(\Gamma =[0,1]\times \{0\}\). We analyze different porous-medium geometries (Tab. 1) and study the influence of the microscale interface roughness, permeability \(\varvec{\textsf{K}}\) and porosity \(\phi\) on the Beavers–Joseph slip coefficient.

To describe a pressure-driven flow, we impose the following boundary conditions for the microscale model (1)–(3):

$$\begin{aligned} \overline{{\varvec{\mathrm v}}} = {\varvec{\mathrm 0}} \quad \text {on }\Gamma _{D},\qquad \overline{{\varvec{\mathrm h}}} = (0,-31.75) \quad \text { on }\,\Gamma _{N}^{\textrm{in}},\qquad \overline{{\varvec{\mathrm h}}} = {\varvec{\mathrm 0}} \quad \text { on }\,\Gamma _{N}^{\textrm{out}}, \end{aligned}$$
(25)

where \(\Gamma _{N}^{\textrm{in}} = \{0\} \times [-0.5,0.5]\), \(\Gamma _{N}^{\textrm{out}}= \{ 1\} \times [-0.5,0.5]\), \(\Gamma _{D} = \partial \Omega {\setminus } (\Gamma _{N}^\textrm{in} \cup \Gamma _{N}^\textrm{out})\). For the coupled macroscale model (4)–(10), we set

$$\begin{array}{cccc} \begin{aligned} \overline{{\varvec{\mathrm v}}}&= {\varvec{\mathrm 0}}{} & {}\quad \text { on }\,\Gamma _{D,{\textrm{ff}}}, \qquad \quad {\overline{v}} = 0{} & {} \quad \text { on }\,\Gamma _{N,{\textrm{pm}}}, \\ \overline{{\varvec{\mathrm h}}}&= (0,- 31.75){} & {} \qquad \text { on }\,\Gamma _{N,{\textrm{ff}}}^{\textrm{in}}, \qquad \quad {\overline{p}} = 31.75{} & {} \quad \text { on }\,\Gamma _{D,{\textrm{pm}}}^\textrm{in}, \\ \overline{{\varvec{\mathrm h}}}&= {\varvec{\mathrm 0}}{} & {} \quad \text { on }\,\Gamma _{N,{\textrm{ff}}}^{\textrm{out}}, \qquad \quad {\overline{p}} = 0{} & {} \quad \text { on }\,\Gamma _{D,{\textrm{pm}}}^\textrm{out}, \end{aligned} \end{array}$$
(26)

where \(\Gamma _{D,{\textrm{ff}}}= \Gamma _{D}\cap \partial \Omega _{\textrm{ff}}\), \(\Gamma _{N,{\textrm{pm}}}= \Gamma _{D}\cap \partial \Omega _{\textrm{pm}}\), \(\Gamma _{N, {\textrm{ff}}}^\mathrm {in/out}= \Gamma _{N}^\mathrm {in/out} \cap \partial \Omega _{\textrm{ff}}\), and \(\Gamma _{D,{\textrm{pm}}}^\mathrm {in/out}= \Gamma _{N}^\mathrm {in/out} \cap \partial \Omega _{\textrm{pm}}\). Boundary conditions (25) and (26) lead to unidirectional flow parallel to the fluid–porous interface, where the normal component of velocity is zero and the pressure field is linear. For this flow problem, the macroscale pressure and the normal velocity component not depend on the choice of the Beavers–Joseph parameter. Therefore, we provide velocity profiles only for the tangential velocity component in the middle of the domain at \(x_1=0.5\).

Below we study different pore geometries. To compute the permeability for geometry \(G_1\), the unit cell problems (17) are solved using an adaptive mesh, where the fluid part \(Y_{\text {f}}\) is partitioned into approximately 35000 elements. For solving the pore-scale problem (1)–(3) in the entire flow domain \(\Omega _{\textrm{mi}}\), an adaptive mesh with approximately 330 000 elements is used. We set the upper boundary \(M=10\) for the Beavers–Joseph slip coefficient in the proposed two-level algorithm (Algorithm 2) in order to consider a broader range of values in comparison with Beavers and Joseph (1967).

4.1.1 Geometry \(G_1\)

In this case, the porous medium is constructed by \(20 \times 10\) periodically distributed circular solid inclusions presented in Tab. 1, which are arranged in line (Fig. 3, left). This yields the characteristic pore size \(\ell =1/20\) and the radius of inclusions is \(r=0.25\ell\). The described pore-scale geometrical configuration leads to a highly porous (\(\phi = 0.804\)) isotropic medium with the permeability tensor \(\varvec{\textsf{K}} = {\text {diag}}\{k_{11}, k_{22}\}\), \(k_{11}=k_{22}\) (Tab. 1). The values of the boundary layer constants appearing in the generalized coupling conditions (11)–(13) are also provided in Tab. 1.

Fig. 3
figure 3

Microscale velocity field (left) and tangential velocity profiles (right) for the pressure-driven flow and geometry \(G_1\)

Full size image

The pore-scale velocity field is presented in Fig. 3 (left). The optimal value of the Beavers–Joseph slip coefficient is \(\alpha _{\textrm{BJ}}^\text {opt} = 2.8\). In Fig. 3 (right), we compare two macroscale Stokes/Darcy problems (4)–(10) with the standard value of the Beavers–Joseph parameter \(\alpha _{\textrm{BJ}}=1\) (profile: SD, \(\alpha _{\textrm{BJ}}=1\)) and the optimal value \(\alpha _{\textrm{BJ}}^\text {opt}= 2.8\) (profile: SD, \(\alpha _{\textrm{BJ}}^\text {opt}=2.8\)) against the pore-scale resolved simulations (profile: pore-scale). In addition, we provide simulation results corresponding to the Stokes/Darcy problem (4)–(7) with the generalized interface conditions (11)–(13) (profile: SD, general). The corresponding errors between the pore-scale and the macroscale simulation results are presented in Tab. 2. The macroscale simulation results for the classical interface conditions with \(\alpha _{\textrm{BJ}}^\text {opt}=2.8\) and the generalized conditions agree significantly better to the pore-scale resolved results than the classical conditions with \(\alpha _{\textrm{BJ}}=1\).

The optimal value of the Beavers–Joseph slip coefficient \(\alpha _{\textrm{BJ}}^\text {opt}\) is independent of velocity magnitude. Several variants of \(\overline{{\varvec{\mathrm h}}}\) and \({\overline{p}}\) in the boundary conditions (25) and (26) have been analyzed, all leading to the same optimal parameter \(\alpha _{\textrm{BJ}}^\text {opt}= 2.8\).

4.1.2 Geometry \(G_2\)

In order to obtain a porous medium with moderate porosity \(\phi = 0.4\), we increase the radius of solid inclusions in comparison with geometry \(G_1\). In this case, we consider \(20\times 10\) circular solid inclusions (\(\ell = 1/20\)) with radius \(r=\ell \sqrt{(1-\phi )/\pi } \approx 0.437 \ell\) arranged in line (Fig. 4, left). This geometrical configuration leads again to an isotropic porous medium and the permeability values together with the boundary layer coefficients are presented in Tab. 1. For geometry \(G_2\), the optimal value of the Beavers–Joseph parameter is \(\alpha _{\textrm{BJ}}^\text {opt}=0.5\) (Tab. 2).

Fig. 4
figure 4

Microscale velocity field (left) and tangential velocity profiles (right) for the pressure-driven flow and geometry \(G_2\)

Full size image

In Fig. 4, we provide velocity profiles corresponding to the pore-scale model, the Stokes/Darcy model with the classical interface conditions (8)–(10) for the two values of the Beavers–Joseph parameter \(\alpha _{\textrm{BJ}}=1\) and \(\alpha _{\textrm{BJ}}^\text {opt}=0.5\), and the Stokes/Darcy model with the generalized coupling conditions  (11)–(13). The corresponding errors between the pore-scale resolved and macroscale solutions are presented in Tab. 2. The classical interface conditions with the optimal value \(\alpha _{\textrm{BJ}}^\text {opt}=0.5,\) as well as the generalized conditions, provide more accurate results than the classical conditions with the traditionally used value \(\alpha _{\textrm{BJ}}=1\).

4.1.3 Geometry \(G_3\)

Here, we study the dependency of the Beavers–Joseph parameter on the pore-scale interface roughness. Therefore, we construct a porous medium which has the same permeability and similar porosity as geometry \(G_1\) but different shape of solid grains. We consider a porous medium consisting of \(20 \times 10\) in-line arranged squared solid inclusions (\(\ell = 1/20\)) of length \(a=0.2154\,\ell\) yielding an isotropic medium (Tab. 1).

We provide the pore-scale velocity field and the profiles of the tangential velocity component in the middle of the domain in Fig. 5. The optimal value of the Beavers–Joseph parameter for this configuration is found by Alg. 2 to be \(\alpha _{\textrm{BJ}}^\text {opt}=7.1\) (Tab. 2). Comparing this result with the one for geometry \(G_1\) (\(\alpha _{\textrm{BJ}}^\text {opt} = 2.8\)) having the same permeability and almost the same porosity, but different pore morphology, we conclude that \(\alpha _{\textrm{BJ}}\) depends not only on the effective properties of the medium but also on the pore geometry including microscale interface roughness. The modified interface condition (14) and the Beavers–Joseph condition (10) with the optimal value \(\alpha _{\textrm{BJ}}^\text {opt} = 7.1\) provide significantly better results than the Beavers–Joseph condition with \(\alpha _{\textrm{BJ}}=1\).

Fig. 5
figure 5

Microscale velocity field (left) and tangential velocity profiles (right) for the pressure-driven flow and geometry \(G_3\)

Full size image

4.1.4 Geometry \(G_4\)

In this test case, we construct a porous medium having the same interface roughness and the same porosity as geometry \(G_1\), but different permeability. For this purpose, we consider \(20 \times 10\) circular solid inclusions arranged in a staggered manner (Tab. 1, Fig. 6). This leads to \(\ell =1/10\) and an orthotropic porous medium with \(\varvec{\textsf{K}} = {\text {diag}}\{k_{11}, k_{22}\}\), \(k_{11}\ne k_{22}\). The radius of the circular solid grains is \(r=0.125\, \ell\). The permeability values and the boundary layer constants in the generalized interface conditions (11)–(13) are given in Tab. 1.

For orthotropic porous media, there exist different approaches in the literature, e.g., (Discacciati et al. 2002; Layton et al. 2003; Eggenweiler and Rybak 2020) to compute \(\sqrt{K}\) appearing in the interface condition (10). As already mentioned in Sect. 2.3, the first one is to take \(\sqrt{K} = \sqrt{{\varvec{\tau }} \varvec{\mathsf {\cdot }}\varvec{\textsf{K}} {\varvec{\tau }}}\), which leads to \(\sqrt{K} = \sqrt{k_{11}}\) for our setting. The second interpretation is \(\sqrt{K} = \sqrt{{\text {tr}}(\varvec{\textsf{K}}) /d }\), where d is the number of space dimensions. In our case, we obtain \(\sqrt{K} = \sqrt{{(k_{11}+k_{22})}/{2}}\). The optimal values of the slip coefficient \(\alpha _{\textrm{BJ}}\) for these two cases are presented in Tab. 2. Note that, taking the second approach for the computation of \(\sqrt{K},\) we obtain the same value \(\alpha _{\textrm{BJ}}^\text {opt}= 2.8\) as for geometry \(G_1\), where the microscale surface roughness is exactly the same.

The choice of \(\sqrt{K}\) influences the optimal Beavers–Joseph parameter \(\alpha _{\textrm{BJ}}^\text {opt}\), however, the ratio \(\sqrt{K}/\alpha _{\textrm{BJ}}^\text {opt}\) is almost the same for both interpretations of \(\sqrt{K}\) with the corresponding value \(\alpha _{\textrm{BJ}}^\text {opt}\). Therefore, in Fig. 6 (right), we present only the macroscale simulation results for the first approach \(\sqrt{K} = \sqrt{k_{11}}\). The microscale velocity field corresponding to geometry \(G_4\) is provided in Fig. 6 (left). We observe that the macroscale simulation results for the classical interface conditions with \(\alpha _{\textrm{BJ}}^\text {opt}=3.0\) as well as for the generalized conditions agree very well to the pore-scale results, however, this is not the case for the classical conditions with \(\alpha _{\textrm{BJ}}=1\).

Fig. 6
figure 6

Microscale velocity field (left) and tangential velocity profiles (right) for the pressure-driven flow and geometry \(G_4\)

Full size image

4.1.5 Geometries \(G_5\) and \(G_6\)

Now we consider two anisotropic porous media with full permeability tensors \(\varvec{\textsf{K}}\). The porous media are composed by \(20\times 10\) elliptical solid inclusions arranged in line (\(\ell =1/20\)) tilted to the right (Tab. 1, geometry \(G_5\)) and tilted to the left (Tab. 1, geometry \(G_6\)). The semi-axes are \(a=0.4 \ell\) and \(b=0.2 \ell\) and the ellipses are rotated clockwise and counterclockwise by \(45^\circ\), respectively. Geometries \(G_5\) and \(G_6\) lead to different permeability tensors (Tab. 1), however, the value of \(\sqrt{K} = \sqrt{{\varvec{\tau }} \varvec{\mathsf {\cdot }}\varvec{\textsf{K}} {\varvec{\tau }}}\) in the Beavers–Joseph condition (10) is the same for our setting. Moreover, both geometries provide the same interface roughness.

Fig. 7
figure 7

Microscale velocity field (left) and tangential velocity profiles (right) for the pressure-driven flow and geometry \(G_5\)

Full size image

We obtain the same optimal value of the Beavers–Joseph slip coefficient \(\alpha _{\textrm{BJ}}^\text {opt}=2.0\) (Tab. 2) for both porous-medium geometrical configurations \(G_5\) and \(G_6\). Therefore, we only present the simulation results for geometry \(G_5\) (Fig. 7). As for geometries \(G_1\) to \(G_4\), we observe that the classical conditions with the optimal Beavers–Joseph parameter \(\alpha _{\textrm{BJ}}^\text {opt}=2.0\) or the generalized interface conditions lead to a more accurate Stokes–Darcy model than the classical coupling concept with the typically used value \(\alpha _{\textrm{BJ}}=1\).

In Tab. 2, we summarize the optimal Beavers–Joseph parameters \(\alpha _{\textrm{BJ}}^\text {opt}\) for all considered geometrical configurations in case of pressure-driven flow. We also present the relative errors (24) for the tangential velocity \(v_1\) in the middle of the domain at \(x_1=0.5\) for the classical interface conditions with two Beavers–Joseph parameters \(\alpha _{\textrm{BJ}}= \alpha _{\textrm{BJ}}^\text {opt}\) (error \(\epsilon _{v_1,\, 0.5, \, \alpha _{\textrm{BJ}}^\text {opt}}\)) and \(\alpha _{\textrm{BJ}}=1\) (error \(\epsilon _{v_1,\, 0.5, \, 1}\)) and for the generalized conditions (error \(\epsilon _{v_1,\, 0.5}\)).

Table 2 Optimal Beavers–Joseph parameters and relative errors for the pressure-driven flow and different porous-medium geometries
Full size table

To summarize, the difference between the pore-scale resolved and the macroscale simulations for parallel flows to the interface can be significantly reduced (Tab. 2) by choosing the correct value of the Beavers–Joseph slip coefficient in the classical interface conditions (8)–(10). Alternatively, one can use the generalized conditions (11), (12) and (14).

4.2 General Filtration Problem

In this section, we study a general filtration problem, where the flow is arbitrary to the porous bed. We show that the Beavers–Joseph condition (10) is not applicable for multidimensional flows to the fluid–porous interface and, therefore, the modified form (14) with the corresponding effective coefficients should be applied. We consider the free-flow region \(\Omega _{\textrm{ff}}=[0,1]\times [0,0.5]\), the interface \(\Gamma =[0,1]\times \{0\}\) and the porous medium \(\Omega _{\textrm{pm}}=[0,1]\times [-0.5,0]\), which includes \(20 \times 10\) circular solid grains (\(\ell =1/20\)) with radius \(r=0.25\ell\) (geometry \(G_1\)). In this case, the interface conditions (8)–(9) and (11)–(12) coincide, since \(N_s^{\textrm{bl}}=0\) (Tab. 1). Therefore, the classical set of coupling conditions and the generalized one differ in the condition on the tangential velocity component: the Beavers–Joseph condition (10) and its generalization (13), which can be written in form of (14).

For the pore-scale problem (1)–(3), we impose the following boundary conditions to obtain arbitrary flow to the porous layer (Fig. 8, left):

$$\begin{aligned} \begin{aligned} \overline{{\varvec{\mathrm v}}}&= (0,-0.7{\text {sin}}(\pi x_1))\, \text { on }\Gamma _{D}^\text {in}\,, \qquad \ \quad \overline{{\varvec{\mathrm v}}} = {\varvec{\mathrm 0}}\, \text { on }\Gamma _{D}^\text {wall} \,, \quad \ \\ \overline{{\varvec{\mathrm h}}} \varvec{\mathsf {\cdot }}{\varvec{\mathrm n}}&= 0\, \text { on } \Gamma ^\text {out}\,, \qquad \overline{{\varvec{\mathrm v}}} \varvec{\mathsf {\cdot }}{\varvec{\tau }} = 0\, \text { on } \Gamma ^\text {out}\,, \end{aligned} \end{aligned}$$
(27)

where \(\Gamma _{D}^\text {in}= [0,1] \times \{0.5\}\), \(\Gamma ^\text {out} = \left( \{0\} \times [0,0.5] \right) \cup \left( \{1\} \times [0,0.225] \right)\), and \(\Gamma _{D}^\text {wall} = \partial \Omega {\setminus } \left( \Gamma _{D}^\text {in} \cup \Gamma ^\text {out}\right)\).

The macroscale Stokes/Darcy model (4)–(7) with the classical interface conditions (8)–(10) or the generalized conditions (11)–(13) is complemented by the following boundary conditions

$$\begin{aligned} \begin{aligned} \overline{{\varvec{\mathrm v}}}&= \left(0,-0.7{\text {sin}}(\pi x_1)\right) \,\text { on }\Gamma _{D}^\text {in} \,, \qquad \ \overline{{\varvec{\mathrm v}}} = {\varvec{\mathrm 0}}\, \text { on }\Gamma _{D,{\textrm{ff}}}^\text {wall} \,, \\ \overline{{\varvec{\mathrm h}}} \varvec{\mathsf {\cdot }}{\varvec{\mathrm n}}&= 0, \quad \overline{{\varvec{\mathrm v}}} \varvec{\mathsf {\cdot }}{\varvec{\tau }} = 0 \qquad \, \text { on } \Gamma ^\text {out}\,, \qquad \ {\overline{v}} = 0 \,\text { on }\Gamma _{N,{\textrm{pm}}}^\text {wall} \,, \end{aligned} \end{aligned}$$
(28)

where \(\Gamma _{D,{\textrm{ff}}}^\text {wall} = \Gamma _{D}^\text {wall} \cap \partial \Omega _{\textrm{ff}}\) and \(\Gamma _{N,{\textrm{pm}}}^\text {wall} = \Gamma _{D}^\text {wall} \cap \partial \Omega _{\textrm{pm}}\).

In Fig. 8 (left), we present the pore-scale velocity field for the general filtration problem with arbitrary flow direction to the porous layer. First, we try to determine an optimal value of the Beavers–Joseph parameter in condition (10) in the same way as in Sect. 4.1. Since the flow is non-parallel to the fluid–porous interface in this case, we consider a broader range for the Beavers–Joseph parameter setting the upper boundary \(M=100\) in Algorithm 2.

In Fig. 8 (right), we provide the pore-scale resolved and macroscale velocity profiles for the tangential component at \(x_1=0.7\), where the flow direction is arbitrary to the interface. The macroscale tangential velocity profile computed using the typical value \(\alpha _{\textrm{BJ}}=1\) does not fit to the pore-scale resolved velocity. We noticed that the difference between the macroscale and the pore-scale simulation results for the tangential velocity \(\epsilon _{v_1, \, 0.7, \, \alpha _{\textrm{BJ}}}\) given by (24) is getting smaller for bigger values of the Beavers–Joseph slip coefficient \(\alpha _{\textrm{BJ}}\), however, these improvements are minor for \(\alpha _{\textrm{BJ}}\in [20,100]\). Therefore, we take \(\alpha _{\textrm{BJ}}^\text {opt}=20\) at \(x_1 = 0.7\) for the tangential velocity. As can be seen from Fig. 8 (right), the macroscale profile for the Beavers–Joseph slip coefficient \(\alpha _{\textrm{BJ}}=20\) fits much better to the pore-scale one than for the traditionally used value \(\alpha _{\textrm{BJ}}=1\).

Since the flow is arbitrary to the interface, we consider four different cross sections and provide the optimal values \(\alpha _{\textrm{BJ}}^\text {opt}\) and the errors (24) for both velocity components \(v_1\) and \(v_2\) in Tab. 3. The microscale surface roughness, porosity and permeability reflected in the Beavers–Joseph parameter \(\alpha _{\textrm{BJ}}\) do not change along the fluid–porous interface for the considered porous media, therefore, the slip coefficient \(\alpha _{\textrm{BJ}}\) must be constant along the interface. However, we observe that the optimal value of the Beavers–Joseph parameter is different for every cross section and each velocity component. This fact indicates the unsuitability of the Beavers–Joseph condition (10) for arbitrary flows.

Fig. 8
figure 8

Microscale velocity field (left) and tangential velocity profiles (right) for the general filtration problem with geometry \(G_1\) and \(\ell =1/20\)

Full size image
Table 3 Optimal Beavers–Joseph parameters and relative errors for the general filtration problem and geometry \(G_1\)
Full size table

In contrast, when the generalized condition (13) or, equivalently, (14) is used instead of condition (10), the pore-scale and the macroscale simulations fit very well (Fig. 8 and Tab. 3). Notice that the errors for the generalized coupling conditions are of the same order as for the classical conditions with the optimal value of the Beavers–Joseph parameter determined at every cross section (Tab. 3). In case of the generalized conditions, all effective coefficients are computed only once. This is a huge advantage in terms of the computational effort. Therefore, the generalized coupling condition on the tangential velocity should be used. We advice the researchers having their codes with the Beavers–Joseph condition to update only two parameters appearing in the generalized condition in form of (14).

5 Conclusions

In this paper, we analyzed the applicability of the Beavers–Joseph interface condition (10) to different flow problems, developed an efficient two-level algorithm to determine the slip coefficient \(\alpha _{\textrm{BJ}}\) and proposed a modification of the Beavers–Joseph interface condition, which is suitable for arbitrary flow directions. We found out that (i) the Beavers–Joseph parameter is in general not constant along the fluid–porous interface and thus the Beavers–Joseph condition is unsuitable for arbitrary flow directions; (ii) the typically used value \(\alpha _{\textrm{BJ}}=1\) is not correct for many flow problems; (iii) the optimal value of the Beavers–Joseph slip coefficient can be found only for unidirectional flows; (iv) the slip coefficient \(\alpha _{\textrm{BJ}}\) depends on the pore-scale geometrical information including the microscale surface roughness; (v) the generalized interface condition (14) should be used instead of the Beavers–Joseph condition (10) for arbitrary flows. These conclusions are made based on the detailed study of two flow problems: the pressure-driven flow parallel to the fluid–porous interface and the general filtration problem, where the fluid flow is arbitrary to the porous bed.

We analyzed different types of porous media: isotropic, orthotropic, and anisotropic. To find the optimal value of the Beavers–Joseph parameter \(\alpha _{\textrm{BJ}},\) we propose an efficient two-level numerical algorithm based on minimizing the difference between the macroscale solution of the Stokes/Darcy problem and the pore-scale resolved solution. In order to study the dependency of the optimal Beavers–Joseph parameter on the pore geometry, we constructed porous media having the same effective properties (porosity, permeability) but different pore morphologies (circular and squared solid grains) and considered different configurations of the porous bed having the same microscale surface roughness.

For parallel flows to the porous medium, the Beavers–Joseph slip coefficient can be determined using pore-scale information that leads to accurate simulation results for the corresponding coupled Stokes/Darcy problem. However, it is not possible to find an optimal value \(\alpha _{\textrm{BJ}}\) for general filtration problems with arbitrary flow directions to the fluid–porous interface. Therefore, the Beavers–Joseph condition is not applicable for such flow problems. In the case of arbitrary flow directions, alternative sets of interface conditions should be used, e.g., those proposed in Angot et al. (2017); Eggenweiler and Rybak (2021); Sudhakar et al. (2021). We reformulated the generalized interface condition on the tangential velocity proposed in Eggenweiler and Rybak (2021) in the form (14) of the Beavers–Joseph condition (10). This modification does not contain any unknown parameters and it is suitable for arbitrary flows in the Stokes/Darcy systems.