A framework for computing effective boundary conditions at the interface between free fluid and a porous medium

Interfacial boundary conditions determined from empirical or ad hoc models remain the standard approach to model fluid flows over porous media, even in situations where the topology of the porous medium is known. We propose a non-empirical and accurate method to compute the effective boundary conditions at the interface between a porous surface and an overlying flow. Using a multiscale expansion (homogenization) approach, we derive a tensorial generalized version of the empirical condition suggested by Beavers & Joseph (J. Fluid Mech., vol. 30 (01), 1967, pp. 197–207). The components of the tensors determining the effective slip velocity at the interface are obtained by solving a set of Stokes equations in a small computational domain near the interface containing both free flow and porous medium. Using the lid-driven cavity flow with a porous bed, we demonstrate that the derived boundary condition is accurate and robust by comparing an effective model to direct numerical simulations. Finally, we provide an open source code that solves the microscale problems and computes the velocity boundary condition without free parameters over any porous bed.

Interfacial boundary conditions determined from empirical or ad hoc models remain the standard approach to model fluid flows over porous media, even in situations where the topology of the porous medium is known. We propose a non-empirical and accurate method to compute the effective boundary conditions at the interface between a porous surface and an overlying flow. Using a multiscale expansion (homogenization) approach, we derive a tensorial generalized version of the empirical condition suggested by Beavers & Joseph (J. Fluid Mech., vol. 30 (01), 1967, pp. 197-207). The components of the tensors determining the effective slip velocity at the interface are obtained by solving a set of Stokes equations in a small computational domain near the interface containing both free flow and porous medium. Using the lid-driven cavity flow with a porous bed, we demonstrate that the derived boundary condition is accurate and robust by comparing an effective model to direct numerical simulations. Finally, we provide an open source code that solves the microscale problems and computes the velocity boundary condition without free parameters over any porous bed.

Introduction
Surfaces found in nature are generally non-smooth with complex hierarchical structural features (Liu & Jiang 2011). The purposes of these surfaces vary greatly, ranging from camouflage and insulation to less obvious functions, such as passively interacting with surrounding fluid to reduce drag or noise (Abdulbari et al. 2013). These functions manifest as effective macroscale properties -for example, permeability, elasticity, slip and optical transparency -while their origin is the small-scale features of the surface. Therefore, to understand the hydrodynamic function of such complex surfaces, a systematic multiscale approach is required. In a bottom-up strategy, the microscale fluid-structure physics of the coating material 868 U. Lācis and S. Bagheri scalar parameters. This approach requires the support of empirical data (Zampogna & Bottaro 2016) or extensive computations to cover a large interval of parameters (Rosti et al. 2015).
In this work, we provide practitioners the framework to compute accurate interfacial velocity boundary conditions, instead of empirically determining them. We derive the interface boundary condition for slip velocity using homogenization and present the relevant Stokes equations to be solved in a microscale interface unit cell. Our main contribution is to provide a set of simple and numerically feasible microscale problems, which once solved, allows for a robust non-empirical effective interface condition. Our interface condition can be considered as a generalized version of the BJ condition, since it depends on interface permeability tensor and on the interface velocity strain rate tensor.
This paper is organised as follows. In § 2, using the lid-driven cavity with a porous bed as an example, we compare the velocity field computed from a direct numerical simulation to the field obtained by solving the homogenized (averaged) equations with the interface condition that is proposed in this paper. After that, in § § 3-5, we derive the interface boundary condition. More specifically, in § 3 we decompose the physical domain into a porous part and free-fluid part and define an interface between the two domains. We then introduce the equations governing the microscale fluid flow in each part as well as their coupling through continuity of velocity and stress at the interface. In § 4, we use multiscale expansion to derive the relevant Stokes equations to be solved in a microscale interface unit cell in order determine the effective boundary condition. In § 5 we derive the interface conditions by employing a homogenization (averaging) technique and relate the obtained results back to the example presented in § 2. Finally, in § 6, we conclude the paper.

Direct numerical simulations versus continuum model
The purpose of this section is to compare two approaches to describe the flow in a lid-driven cavity with a homogenous bed of solid cylinders. In the first approach, represented in figure 1(a), we solve the Stokes equations over all spatial scales. This is possible for simplified geometries such as this one, but clearly for more complex (possibly three-dimensional) and much denser porous beds it is not feasible to solve Stokes equations with complete microscale resolution. In the second approach represented in figure 1(c), we reduce the degrees of freedom of the flow in the porous bed by homogenization.
The cavity has a length H and a depth of (H + d), where the porous bed is confined to −d < y < 0 and −H/2 < x < H/2. Coordinate y = 0 corresponds to the tangent plane of the top row of cylinders and the depth of the porous bed is d ≈ H/2. The top wall of the cavity is driven by a constant streamwise velocity U w , which is sufficiently slow to render fluid inertia negligible inside the cavity. Figure 1(b) shows the geometry of the porous bed, which consists of a lattice of cylinders with diameter D = 2r and with the spacing l. For the particular example discussed in this section, the microscale length is l/H = 0.1 and the cylinder volume fraction is φ = πr 2 /l 2 = 0.02.

Direct numerical simulations
The two-dimensional Stokes equations are solved with a no-slip condition imposed at the cylinder surfaces as well as on the vertical and bottom walls of the cavity. The Left frame (a) shows the lid-driven cavity domain with a porous bed. The top wall is driven with velocity U w . Centre frame (b) shows a magnified view of the porous bed. Cylinders are spaced apart by distance l and Γ c defines the boundary of the cylinders. Right frame (c) shows a two-domain description of the same cavity problem in a homogenized setting. The parameters defining the problem are the volume fraction φ = πr 2 /l 2 = 0.02, scale separation l/H = 0.1 and porous-bed depth of d ≈ 0.5H. equations are given by, where µ is fluid viscosity, Γ n denotes the bottom and side boundaries of the cavity and Γ c denotes the boundary of the cylinders in the porous bed. The computations are performed with FreeFEM++ (Hecht 2012), using a triangular mesh and a Taylor-Hood finite element space (P2+P1) for velocity and pressure. We set mesh spacing to s 1 = 0.125l at the outer boundaries (cavity walls) of the domain, and s 2 = 0.050l at the surface of the cylinders. We have carried out a simulation of the same configuration with half the mesh spacing ( s 1 = 0.063l and s 2 = 0.025l), and observed that the slip velocity changes by 0.6 %. In figures 2 and 3, we present the velocity profiles obtained from DNS (direct numerical simulations) with solid black lines. Figure 2 shows the streamwise and wallnormal velocity profiles for the fixed streamwise position x/H = −0.1. The insets show the detailed microscale fluctuations of the velocities near the interface and the rapid transition to the macroscale velocity in the free-fluid region. This transition occurs in a thin layer near the top row of cylinders. Figure 3 shows the velocity along x at a virtual free-fluid-porous interface placed at y s /l = 0.1. If y s would have been an interface with a rigid wall, these graphs would show the no-slip condition, i.e. zero velocity for both wall-normal and streamwise velocity components. However, at the interface with a porous medium, there is a slip velocity u s and a penetration velocity v p . The underlying structure of porous medium manifests as microscale oscillations in slip and penetration velocities. Apart from the microscale oscillations, we can observe that both velocity components exhibit 870 U. Lācis and S. Bagheri  macroscale variations. The negative slip velocity, which is induced by the spanwise vortex above the interface, has a maximum value u s /U w = −0.0121 at the centre of the cavity. Although the slip velocity is small compared to the bulk flow, it may have a significant physical effect on the characteristics of the overlying fluid.
For example Rosti et al. (2015) recently showed that slip velocities below 3 % had a significant effect on the flow statistics in a turbulent channel. Additionally, Carotenuto & Minale (2013) showed that precise predictions of slip velocity can be essential to get accurate viscosity measurements from rheology tests. The penetration velocity shows a sinusoidal behaviour; for x < 0, there is a net mass transport from the pore region to the free-fluid region, whereas for x > 0 the net mass flow is in the opposite direction. Similarly, the net momentum transport into the porous region is in opposite directions whether x > 0 or x < 0. The values of u s and v p , which are essential to capture the momentum and mass transport across the interface, depend both on the flow in the pores and on the microscale geometry of the pores. In the next section, we introduce a fully non-empirical method to compute u s and v p with an error of O(l/H) without resorting to DNS of the full domain.

Simulation of homogenized equations
We start by replacing the full DNS domain with two rectangular domains, where the free-fluid region Ω f and porous region Ω p are separated by the interface Γ , as shown in figure 1(c). In the free-fluid part, we do not employ homogenization and therefore the flow is governed by the Stokes equations whereû andp are flow and pressure fields in Ω f , respectively. Dirichlet conditions are imposed forû on the vertical side walls and the top wall of the cavity (as for the full DNS). The boundary condition at the interface Γ in contact with the porous region isû = (û s ,v p ), (2.7) where the slip velocity and the penetration velocity depend on the flow in the porous region. In the porous part Ω p , the flow is governed by the well-known Darcy's laŵ and mass conservation where K itr is the interior permeability tensor. It is convenient to combine mass conservation with Darcy's law to arrive at a single equation for the pore pressure, which -assuming that permeability tensor K itr is constant over space and isotropicreads p = 0. (2.10) We complement this equation with homogenous Neumann conditions on the side walls and the bottom wall, which corresponds to zero transpiration. At the interface Γ , we impose a Dirichlet condition with the pressure obtained from (2.5)-(2.6). This continuous pressure condition is valid up to O(l/H) under the theoretical assumptions employed in this paper, which will be discussed in following sections. Since the particular porous bed we are investigating is isotropic (see figure 1a,b), the continuity of pressure at the interface is also implied by works of Marciniak-Czochra & Mikelić (2012) and Carraro et al. (2013).

U. Lācis and S. Bagheri
Returning to the velocity boundary conditions (2.7), that are required for solving Stokes system (2.5)-(2.6), we simply state the conditions of order O(l/H) that will be derived in the next sections (with the final result in (5.11)). The penetration velocity componentv p is given byv where K cyl is the isotropic permeability of the porous medium consisting of a regular array of circular cylinders. Note that, although pressure is continuous at Γ in our case, the pressure gradient is not necessarily continuous; ∂ yp − in (2.11) denotes the pressure gradient when approaching the interface from the porous bed. The condition for v p can also be obtained from mass conservation for a thin rectangular control volume around Γ with periodic streamwise velocity on the vertical sides. The condition for slip velocityû s isû This expression is similar to the condition obtained empirically by Beavers and Joseph, except that K s is the interface permeability (e.g. K s = K cyl ), related to a semi-permeable transition layer between the porous medium and the free fluid.
Another difference with the BJ condition is that the strain term ∂ xv is included in addition to ∂ yû . It has been argued that the term ∂ xv should be present for curved boundaries (Jones 1973), but to the authors' knowledge, it has not been derived earlier for flat interfaces, although it has be conjectured to exist by Nield (2009). The constant L s is related to the slip length in the Navier boundary condition. The constants appearing in boundary conditions (2.11)-(2.12) are provided by microscale simulations in interface cells, described in following sections. In order to provide an overview of the applicability of the derived boundary conditions, we summarize here briefly the practical limits that will be determined both theoretically and numerically in the remaining part of this paper. These limits are; (i) moderate scale separation l/H 0.1; (ii) restriction on the Reynolds number based on the seepage velocity (Re d 1); and (iii) restriction on the Reynolds number based on the lid velocity U w (Re f −1 ). The corresponding Reynolds numbers are defined later. We solve the set of equations (2.5)-(2.12) using FreeFEM++ with mesh spacing s = 0.125l. Figure 2 (blue curve with circular symbols) compares the obtained velocity profiles over a vertical slice to the DNS results, where one can observe an excellent agreement between the two. We note that the effective macroscale behaviour is captured, while underlying oscillations arising from the small-scale characteristics of porous bed are not modelled. The consequence of using an O(l/H) accurate model in the interior (Darcy's law) is that the diffusion process from the free flow to the pore flow (which defines the transition layer of height ∼ l) is not captured. However, from the perspective of the free fluid, the macroscopic effect of the porous bed is essentially the same using fully resolved DNS and the continuum model. In figure 3(a) (blue curve with circular symbols), we also observe a good agreement between DNS and the continuum model for slip velocity over a horizontal slice, despite that the latter approach does not resolve the microscale dynamics between and around the cylinders. The model is able to predict the maximum slip velocity at the centre of the cavityû s /U w = −0.0119. Figure 3(b) compares the predicted penetration velocity and DNS results. There one can observe that, although microscale oscillations dominate the DNS results, the macroscale sinusoidal behaviour is correctly captured by the model. To close this section, we want to point out that the effective problem is computationally much cheaper than DNS. The number of degrees of freedom used for the DNS in § 2.1 for the region below the interface is approximately 2.0 × 10 5 , whereas for the two-domain approach in § 2.2 it is 1.4 × 10 4 . This difference arises from coarser mesh in the porous region, as well as the reduced number of variables (the model equation defines pressure only). In more complex and three-dimensional cases the difference can be significantly larger, and only averaged models might be computationally feasible to solve numerically.

Governing equations and flow decomposition
While the derived effective boundary condition in this paper has been applied to the steady cavity flow over regular array of cylinders, the boundary condition is more general and can be applied to more complex three-dimensional flows, as schematically shown in figure 4(a). We relax the assumptions of a steady and non-inertial flow in the cavity, and start from incompressible Navier-Stokes equations. The length H now corresponds to an appropriate macroscopic length scale of the flow, satisfying = l/H 1.

Dimensionless Navier-Stokes equations
The free-fluid region and the porous region are characterized by different spatial and temporal scales. We choose to non-dimensionalize the Navier-Stokes equations using the characteristic scales of the porous medium. In appendix A, the reader can find detailed analysis of these scales. The characteristic velocity of the flow in the porous region U d is where P is the characteristic macroscopic global pressure, µ is the fluid viscosity and H and l are the macroscopic and microscopic length scales, respectively.

U. Lācis and S. Bagheri
Consequently, we use the following relationships between dimensional (denoted with 'tilde') and dimensionless variables Here, time is non-dimensionalized with the convection time scale at the microscale. In order to simplify the notation, we use x 1 and x, and x 2 and y interchangeably. We may now write the Navier-Stokes equations in the following dimensionless form where Re d = ρ f U d H/µ is the Darcy Reynolds number. The different order of the scale separation parameter in front of the terms provides an estimate of the relative magnitude of the terms within the porous region; pressure force plays a dominant role and the inertial force is a higher-order effect. This conclusion holds if Re d 1, which is an assumption for the presented approach. Technically, equations (3.3)-(3.4) hold also in the free-fluid region, although the relative magnitude between the terms is not anymore characterized by the scale separation parameter; the scaling (3.2), except for pressure, is not suited for the free flow.
3.2. Decomposition of the flow field We continue by choosing an interface Γ at a vertical coordinate x 2 = y s , which divides the fluid domain into a free-fluid region and a porous region. The accuracy of the final interface condition does not depend on y s (up to certain limits) as proven by Marciniak-Czochra & Mikelić (2012). We confirm this statement numerically in the § 5.3.
We separate the flow above the interface (domain Ω f ) into a fast flow (U i , P) and a perturbation (u + i , p + ), The terms (u + i , p + ) are generated by the porous medium and will -as shown below -be responsible for the induced slip and penetration velocities. The pressure and the velocity below the interface (domain Ω p ), denoted by represent the slow flow and the pressure field in the pores. Table 1 summarizes the introduced quantities in Ω f and Ω p . By inserting the decomposition (3.5) and the quantities (3.6) into (3.3)-(3.4) and grouping the different terms, the equations governing the dynamics of the different quantities are obtained.

Fast flow
The global pressure P and the fast flow U i are governed by the Navier-Stokes equations with no-slip condition at y s , i.e.
Boundary conditions at interface between free fluid and a porous medium Using the introduced non-dimensionalization (3.2) on the variable estimates, we arrive at a priori orders of the fast flow velocity and the global pressure (3.10) The dimensionless fast flow field U i becomes very large for small , whereas the global pressure in Ω f -either externally imposed as in pressure-driven channel flow or induced by the fast flow as in lid-driven cavity -is of order one. Note that the assumption (3.9) is not a universal one, and depends on the bulk Reynolds number; for example, for Stokes flow U f ∼ −2 U d . Appendix A shows that (3.9) is obtained, when Re f = ρ f HU f /µ ∼ −1 . This a priori scale estimate simplifies the multiscale expansion outlined in § 4 and its consequence is of a theoretical nature; in practice, our derived boundary condition predicts very accurately the slip and penetration velocity for a wide range of parameters, as we demonstrate in § 5.3.

Perturbations and slow flow
The perturbations above the interface are governed by In order to solve these equations, one has to know the flow below the interface (u − i ). The pressure and the slow velocity below the interface are governed by where the stress tensors containing u + i , u − i and U i in (3.14) are defined by (B 7)-(B 9) in appendix B. Note that the decomposition of the flow introduced in § 3.2 is exact, since continuity of both velocity (3.12) and total stress (3.14) is imposed at the interior interface Γ . In other words, if one would sum up (3.7), (3.11), (3.13) and the boundary conditions (3.8), (3.12), (3.14), the Navier-Stokes equations (3.4) defined in the full domain would be recovered.
The perturbation terms (u + i , p + ) in Ω f are an effect of the porous medium. Therefore the perturbation velocity is estimated by the seepage velocity, in dimensional setting u + i ∼ U d , and the pressure perturbation by the microscale pressurep + ∼ p, where p is the pressure induced by the seepage velocity, see appendix A. In the current non-dimensional setting, we have Comparing to (3.10), we observe p + P and u + i U i when 1. The slow velocity below the interface can be estimated asũ − i ∼ U d , since it is reasonable to assume that the flow inside the porous medium has the magnitude of the Darcy velocity. The pore pressure, on the other hand, can be estimated asp − ∼ P, because the global macroscale pressure is present also in the porous medium. For dimensionless variables we then have The dimensional estimates and dimensionless orders of the decomposed quantities are summarized in table 1.

Multiscale expansion
We now turn our attention to the multiscale analysis of the flow near the interface and construct an approximate description of it within an interface cell (defined below). To carry out the multiscale expansion, we introduce the macroscale and microscale coordinates respectively. These coordinates are appropriate to describe the macroscopic and microscopic variations and are related to each other by X i = x i . In the new coordinates, there are two derivatives appearing due to the chain rule where () ,i 0 denotes the derivative with respect to X i and () ,i 1 with respect to x i . The fast flow U i and the global pressure P do not depend on microscale coordinate, i.e. U i,j 1 = 0 and P ,i 1 = 0. This is a direct consequence of the definition of fast flow problem (3.7)-(3.8) and is valid for 1. For u ± i and p ± , which depend on both coordinates, we carry out the multiscale expansion as explained by Mei & Vernescu (2010). The perturbation velocity and the pressure above the interface (y y s ) are expanded as (4.4) Boundary conditions at interface between free fluid and a porous medium

877
The pressure expansion starts with the O( ) term, since p + = O( ). Below the interface (y y s ) the slow flow and the pore pressure are expanded as, (4.6) We insert expansions (4.3)-(4.6) into the corresponding equations (3.11)-(3.14), and collect the terms at different orders. In the following subsections, we introduce and solve equation systems appearing at first two orders (O(1) and O( )).

O(1) equation and its analytical solution in an interface cell
Collecting the terms with prefactor 1, we get the following system We observe that the zeroth-order pressure in the porous region p −(0) is independent of the microscale coordinate x i . For our purpose, which is to derive the macroscale effective boundary condition, it is sufficient to solve this equation in an elongated cell near the vicinity of the interface (figure 4b). The size of this cell is chosen in such a way to capture only the microscale behaviour near the interface. Recall that we use x to describe the streamwise coordinate x 1 and y to describe the interface-normal coordinate x 2 . The solution to (4.7)-(4.8) below the interface in Ω cell is constant and equal to global pressure P at the interface, p −(0) = P| y s , y 1 y y s . (4.10) We will use this result in § 5.1 to derive a macroscale pressure condition at the interface.

O( ) equation and its computational solution in an interface cell
Next, we collect the first-order terms with prefactor , which results in the Stokes equations for (u ±(0) i , p ±(1) ). Specifically, above the interface, we have where A 1 denotes the Stokes operator, which contains derivatives with respect to microscale x i , as defined in (B 4). Below the interface, we have We observe that the slow flow u −(0) i is forced by the macroscale gradient of the pore pressure term (p −(0) ,i 0 ). In contrast, above the interface the equation for perturbation u +(0) i is not driven by a lower-order pressure term; the O(1) pressure above the surface is P, which is contained in equations for U i (3.7)-(3.8). The same global macroscale pressure is driving the fast flow and the slow flow, but whereas this pressure is defined 878 U. Lācis and S. Bagheri as P above the interface, it is obtained as the leading-order expansion term below the interface. This is a consequence of the fact that the fast flow is not expanded, since it varies only with macroscale coordinate by definition.
The boundary conditions of O( ) equations at the interface y s are continuity of velocity and jump in stress (4.14) Here, Σ u ± (1) ij = −p ±(1) δ ij + (u ±(0) i,j 1 + u ±(0) j,i 1 ) is the stress tensor of the perturbation velocity and the slow velocity, whereas S ij = U i,j + U j,i is the strain tensor of the fast flow.
The system of equations (4.11)-(4.14) is solved in the interface cell Ω cell shown in figure 4(b). To complete the problem formulation, boundary conditions are needed at the sides of the cell. Due to regularity of the porous structure, we impose periodic conditions on the sides, as shown in figure 4(c). At the bottom of the cell, we impose the interior solution, which is where K itr ij is the classical interior permeability field. For the derivation of this expression and the corresponding microscale problem, the reader is referred to the book by Mei & Vernescu (2010). From the literature (Mikelić & Jäger 2000;Jäger & Mikelić 2009;Marciniak-Czochra & Mikelić 2012;Carraro et al. 2013) it is well known that the interface cell is exposed to a zero-stress condition at the infinity within the method of matched asymptotic expansion, therefore we impose zero-stress condition, i.e. Σ u + (1) ij n j = 0, at y = y 2 . We may consider the Stokes equations (4.11)-(4.12) and the boundary conditions (4.13)-(4.15) as one linear problem with four unknowns (u ±(0) i , p ±(1) ). Due to linearity, we can construct the solutions for the velocity and pressure fields as the superposition of p −(0) ,i 0 (X i ) and S ij (X i ), i.e. u ±(0) (4.16) and The average of K ij is a tensorial effective Darcy permeability and the average of L ijk is related to the tensorial version of the Navier-slip coefficient. Similarly, the tensors A j and B ij are transfer coefficients from the driving pressure gradient and fluid strain, respectively, to the perturbation pressure.

Microscale Stokes problems for Darcy term
By inserting the ansatzes (4.16)-(4.17) into (4.11)-(4.15), it follows that the tensors K ± ij and A ± i satisfy, (4.19) with boundary conditions at the interface y s given by (4.20) Boundary conditions at interface between free fluid and a porous medium At the bottom boundary y = y 1 , we have K − ik = K itr ik (note that fields K − ik and K itr ik are not coinciding, but have the same values at y = y 1 ). The field K ± jk represents the jth velocity component of the kth Stokes problem. Thus to determine every component of K ± ij and A ± i , 3 pairs of Stokes problems have to be solved coupled at the interface through continuity of velocity field and stress. Note that below the interface, the flow is driven by a unit forcing in one direction at a time. Therefore the physical interpretation of K i1 for example, is the flow response to forcing in the horizontal direction below the interface, as shown in figure 5(a,b).
To obtain reliable results, the interface cell needs to extend sufficiently into the free fluid such that variations of K + ij are small and sufficiently into the porous medium such that variations of K − ij are periodic. We have investigated different heights of the interface cell, and have determined that height of 10l (containing 5 cylinders below the interface) is sufficient. We have computed the K ij and L ijk fields for finer meshes and found no numerical artefacts nor noticeable modification of the results. Figure 5 shows K ± ij fields and corresponding plane-averaged profiles near the tip of the solid structure for interface location y s = 0.1.

Microscale Stokes problems for Navier-slip term
By inserting the ansatzes (4.16)-(4.17) into (4.11)-(4.15), the following equations for the tensors L ± ijk and B ± ij are obtained, (4.22) with boundary conditions at y s given by, Stokes problems. The forcing for these equations is at the interface in the form of a stress condition. For example, the L i12 component is the flow response to a unit tangential stress at the interface, whereas L i22 is the response to unit normal stress at the interface. In general, for problems with flat interfaces described in a coordinate system aligned with the interface, only three pairs of problems are forced. Returning to our two-dimensional configuration with a flat interface and aligned coordinate system, L ij1 are unforced problems, leading to trivial solutions for all components. Out of the forced problems L ij2 , the components L 122 and L 222 are zero since the forcing is in a constrained direction, i.e. due to mass conservation, the motion in the vertical direction is zero when the no-slip condition is enforced at the bottom boundary. We are thus left with only one non-trivial problem (L i12 ), for which the flow fields are shown in figure 6, along with corresponding plane-averaged profiles.

Effective interface conditions
This section provides the final forms of the effective boundary conditions by averaging the microscale solutions provided in the previous section. Since the interface cell is located at the boundary between the free-fluid and porous regions, the conditions for the free fluid should be evaluated using the solution above the interface, while conditions for the porous region should be evaluated using the solution below the interface. In particular, to solve for the pore pressure Laplace equation (2.10) below the interface, one needs a boundary condition for the pressure at the interface; this can be obtained by averaging the O(1)-problem (see § 4.1). To solve the velocity of the free flow above the interface, one needs a condition for the velocity at the interface. We recall however that the O(1)-problem given by (4.7)-(4.8) does not contain velocity. Therefore, one has to investigate solution of the O( )-problem (4.16) to determine a boundary condition for the velocity. This is a consequence of the pore velocity viscous term being higher order compared to the pressure gradient term, see prefactors in (3.3).
5.1. Condition for pore pressure Using decomposition (3.5) and the scaling for pressure perturbation (3.15), we can write the pressure field above the interface as p = P + O( ). (5.1) Taking the average of the expression above in a microscale volume of size l 3 -see definition (B 17) -above the interface gives, Here, we have no angle brackets around P, since it is independent of the microscale coordinate. The multiscale expanded pressure field (4.6) below the interface, on the other hand, is where its volume-averaged form is According to solution of the O(1)-problem (4.10), one can state that the pressure field in the whole interface cell is constant and equal to the macroscale pressure P.
Inserting the solution of the O(1)-problem (4.10) in (5.4) and then equating it to expression (5.2), we obtain the following at the interface, For brevity, we denote averaged dimensional quantities with a 'hat' (e.g.p − = P p − ), which gives the pressure interface condition in its final dimensional form, Working with the chosen estimates (see table 1), one obtains pressure continuity up to O( ) for any anisotropic porous bed. We point out that from (4.17), one may formulate a pressure condition valid to O( 2 ). This is however out of the scope for this work. We note that this result is different compared to works by Marciniak-Czochra & Mikelić (2012) and Carraro et al. (2013). This is a direct consequence of the theoretical assumption Re f ∼ −1 , which leads to the O(1)-problem for pressure being trivial.

Velocity boundary condition for free fluid
The solution to the O( )-problem (see § 4.2) is obtained numerically by computing K ij and L ijk . One may then proceed to construct the fully resolved flow field with error O( ) near the interface. First, we write the velocity above the interface as The macroscale term U i does not appear in the expression above, because it is constant in the interface cell and this constant has be zero due to the boundary condition (3.8).
Inserting the above expression into (4.16) gives

U. Lācis and S. Bagheri
We average out the microscale oscillations by forming the volume average above the interface, Here, we have no angle brackets around pressure and velocity gradients since these quantities are independent of the microscale coordinate. Now, by inserting the approximations of the velocity strain (C 3) and the pressure gradient (C 4) -see appendix C -into (5.9), we obtain where for convenience, we have denoted K ij = K + ij and L ijk = L + ijk . The coefficients K ij and L ijk are evaluated over l 3 cube (l 2 in 2D setting) at the top of the finite interface cell. Finally, we have used that u i = U i + u + i = u + i at the interface. In order to return to the boundary conditions (2.11)-(2.12) used for the lid-driven cavity problem in § 2, we revert the boundary condition to dimensional quantities, The coefficients, based on solutions of the interface cell problems, for the cavity flow are K cyl = l 2 K 22 , K s = l 2 K 11 and L s = l L 112 . (5.12) All other components of tensors K ij and L ijk are zero for the porous bed with regular circular cylinders, and therefore there is no velocity shear term for the penetration velocity (2.11). Actual values used in § 2 are K 11 = 0.0312, K 22 = 0.0986 and L 112 = 0.1783. Equation (5.11) is the final expression of the velocity boundary condition for a rigid porous bed. It can be used together with Navier-Stokes equations in any domain of interest, in order to take into account the effects of the porous medium, without resolving the microscale flow within the porous bed. We emphasize that the 'minus' notation for pressure means that the pressure gradient in the boundary condition is the gradient of the pore pressure. In the final subsection, we test the robustness of the derived boundary condition by varying solid volume fraction, scale separation parameter and interface location.

Accuracy of slip prediction and robustness to interface location
We now return to the example of the lid-driven cavity with a porous bed in order to illustrate more quantitatively that the proposed boundary condition yields accurate and robust slip velocity predictions. More specifically, we carry out a parametric study and report predictions of maximum slip velocityû s at the centre of the cavity. In order to do a fair comparison, we surface average the DNS resultsū s = u s S at the interface. We also assess the contributions from two different terms in the derived boundary condition (5.11). Table 2 shows that for the range of parameters considered, the contribution at the interface from the Navier-slip term is always at least an order of magnitude larger than the contribution from the Darcy term. This is a consequence of (5.11), where K s ∼ l 2 and L s ∼ l, and therefore K s L s for fine microstructures. This result is in agreement with previous work. It was first suggested by Saffman (1971) that the Darcy term is of higher order and can be neglected, the same result was y s φû sK /ū sûsL /ū s (×10 −1 )û s /ū sūs /ū s later rigorously proved by Mikelić & Jäger (2000). One may therefore obtain a good approximation with only the Navier-slip term, as first suggested by Saffman (1971) and later rigorously shown by Mikelić & Jäger (2000). Including the Darcy term however yields -consistently -a smaller error, and therefore also a robust velocity boundary condition with respect to the interface location and different pore geometries. Additionally, the Darcy term is the only contribution appearing in the interface-normal direction, which is essential to capture the momentum transfer from and to the porous region.
Note that although the Darcy term is much smaller than the Navier-slip term, both terms appear in the O( ) equation (see § 4.2). This is a consequence of the estimate that perturbation velocity is of the same order as the Darcy velocity, i.e. u ± i ∼ U d and that the fast flow velocity scales as (3.9). An alternative approach would be to assume U f ∼ −2 U d , which would essentially result in three velocity scales to allow for the slip velocity to be faster than the Darcy flow but slower than the free flow. In such an approach -which would more sophisticated than the current one -the Darcy term can appear in higher-order equation than the equation for which the Navier-slip term appears. Nevertheless, one can observe that our direct numerical simulations are in good agreement with simulations of the homogenized model despite that we strictly speaking do not satisfy (3.9) for Stokes flow. Thus our method can be considered as practical parameter-free framework for computing the coefficients of a generalized condition proposed by Beavers and Joseph; a condition which has been employed under a variety of different flow conditions by experimentalists. To sum up, we have theoretically assumed that (i) 1; (ii) Re d 1 and (iii) Re f ∼ −1 . Based on numerical tests in this section for Re f = 0 and up to = 0.1, we have determined the practical limits of the derived interface condition by relaxing conditions (i) and (iii), as summarized in § 2.2.

Conclusions
We have presented a framework to construct a reduced homogenized model of the flow above and through a porous medium consisting of regular solid structures of general shape. The main contribution of the present paper is to provide the foundation and the tools to compute effective boundary conditions completely free from data fitting. The approach that we adopted can be summarized by the following four steps. First, the governing equations are made non-dimensional using scale estimates arising from flow in the porous domain. Second, the governing equations describing the fully resolved flow are separated at a virtual interface, and decomposed above the interface into equations for the fast flow and for perturbations. Third, multiscale expansion according to Mei & Vernescu (2010) is employed on perturbation and pore equations. Finally and after solving O(1)-and O( )-problems, we construct the interface conditions for pore pressure and fluid velocity using volume averages.
This procedure results in macroscopic description of the flow over a porous domain with an error O( ). Specifically, using our a priori scaling estimates, the leading-order conditions are a pressure continuity condition and a generalized tensorial BJ boundary condition. The proposed velocity condition depends on the interface permeability and on the velocity strain, while the BJ condition contains interior permeability and velocity derivative of one component in one direction only. To the authors' knowledge, such a general formulation has not been derived and validated before. Moreover, in order to obtain the constants of the effective boundary conditions, we derive a number of Stokes problems that need be solved numerically in small interface unit cells. Solvers for microscale problems have been released as an open source software (Lācis & Bagheri 2016), along with the solver used for the lid-driven cavity flow.
This work is also among the first to validate non-empirical boundary conditions on two-dimensional flows with DNS, where penetration velocity, slip velocity and pressure conditions have to be predicted to solve the coupled two-domain problem. The present boundary condition has been tested in the lid-driven cavity flow for a range of volume fractions from φ = 0.02 to 0.45, scale separation parameters from = 0.02 to 0.1 and interface locations from y s = 0.1 to y s = 0.5. When the homogenized model results are compared to DNS, the slip velocity predictions have been found to be robust and to give accurate predictions for all investigated parameters. We hope that, with this work and the release of the associated software, we can provide the numerical fluid dynamics community the tools to model flows over existing non-smooth surfaces as well as to design surfaces to modify fluid flow characteristics. gradient at the pore scale and ∇ = () ,j 0 denotes the gradient at the macroscale. This is a classical result at first order in the interior, as derived by Mei & Vernescu (2010) and also used by Gopinath & Mahadevan (2011). The force balance is thus where U d is the characteristic velocity in the porous region and p is characteristic microscopic pressure.
Comparing the first and second terms, we immediately arrive at the first estimate used in the main paper (3.1). We argue in the main paper that the perturbation velocities are caused by the porous medium and therefore could be estimated based on the characteristic velocity in the porous regioñ The pore pressure is associated with the global pressure difference and estimated as p − ∼ P, (A 4) whereas the pressure perturbation above the interface we associate with the microscopic pressure difference p. We argue that the pressure perturbation above is caused directly by the flow in the pores, which results in pressure difference in the pore scale. The estimate we use isp + ∼ p.

(A 5)
For the fast flow, we use an estimateŨ where U f is the characteristic fast flow velocity. The associated pressure we estimate using the same macroscopic pressure difference as everywhere elsẽ The magnitude of the fast flow velocity U f (3.9) can be estimated from momentum balance in the free fluid. The momentum of the fast flow is governed by ρ f (∂ t U + (U · ∇)U) = µ∇ 2 U − ∇P.

886
U. Lācis and S. Bagheri where Re f = ρ f HU f /µ is free-fluid Reynolds number. Now one has to say something about Re f in order to finalize estimate of U f . One possible choice is to assume Re f = O( −1 ), and then we recover assumption (3.9). The estimate of U f is later implicitly used in order to determine the order of shear stress from the free fluid at the interfacẽ where we have assumed that the fast free-flow velocity is obtained over the macroscopic length scale. Finally, estimates (A 3)-(A 7) can be made non-dimensional following section (3.2). After using the momentum balance presented here (A 2) and assumption (3.9), one arrives at the dimensionless orders presented in rightmost column of table 1. Additionally, making the shear-stress estimate (A 12) non-dimensional, one obtains U i,j = O(1), which together with prefactor in the non-dimensional stress (B 9) leads to free-fluid shear appearance in O( )-problem (4.14). and stress tensors for the Navier-slip interface problem are Σ L + ij = −B + kl δ ij + (L + ikl,j 1 + L + jkl,i 1 ), (B 14) The velocity strain tensor is S ij = U i,j + U j,i . f (x − X) dx dy, (B 18) where the plane is oriented parallel to the interface. Here, X is the location of the centre of the averaging surface. In the two-dimensional case, as in the lid-driven cavity example reported in the paper, this integral reduces to a line integral. For convenience, the X argument is omitted in the main paper.