The asymptotic analysis of a Darcy-Stokes system coupled through a curved interface

The asymptotic analysis of a Darcy-Stokes system modeling the fluid exchange between a narrow channel (Stokes) and a porous medium (Darcy) coupled through a $ C^{2} $ curved interface, is presented. The channel is a cylindrical domain between the interface ($ \Gamma $) and a parallel translation of it ($ \Gamma + \epsilon \, \boldsymbol{\widehat{e}}_{N} $). The introduction of a change variable to fix the domain's geometry and the introduction of two systems of coordinates: the Cartesian and a local one (consistent with the geometry of the surface), permit to find a Darcy-Brinkman lower dimensional coupled system as the limiting form, when the width of the channel tends to zero ($ \epsilon \rightarrow 0 $).


Introduction
In this paper we continue the work presented in [14], extending the result to a more general scenario. That is, we find the limiting form of a Darcy-Stokes (see (26) ) coupled system, within a saturated domain Ω in R N , consisting in three parts: a porous medium Ω 1 (Darcy flow), a narrow channel Ω 2 whose width is of order (Stokes flow) and a coupling interface Γ = ∂Ω 1 ∩ ∂Ω 2 , see Figure 1 (a). In contrast with the system studied in [14], where the interface is flat, here the analysis is extended to curved interfaces. It will be seen that the limit is a fully-coupled system consisting of Darcy flow in the porous medium Ω 1 and a Brinkman-type flow on the part Γ of its boundary which now takes the form of a N − 1 dimensional manifold.
The central motivation in looking for the limiting problem of our Darcy-Stokes system is to attain a new model free of the singularities present in (26). These are the narrowness of the channel O( ) and the high velocity of the fluid in the channel O( ); both (geometry and velocity) with respect to the porous medium. Both singularities have substantial negative computational impact at the time of implementing the system, such as numerical instability and poor quality of the solutions. Moreover, when considering the case of curved interfaces, the geometry of the surface intensifies these effects, making even more relevant the search for an approximate singularity-free system as it is done here.
The relevance of the Darcy-Stokes system itself, as well as its limiting form (a Darcy-Brinkman system) is confirmed by the numerous achievements reported in the literature: see [4], [2], [6] for the analytical approach, [3], [5], [9], [13] for the numerical analysis point of view, see [11], [21] for numerical experimental coupling and [12] for a broad perspective and references. Moreover, the modeling and scaling of the problem have already been extensively justified in [14], hence, this work is focused on addressing (rigorously) the interface geometry impact  in the asymptotic analysis of the problem. It is important to consider the curvature of interfaces in the problem, rather than limiting the analysis to flat or periodic interfaces, because the fissures in a natural bedrock (where this phenomenon takes place) have wild geometry. In [7], [8] the analysis is made using homogenisation techniques for periodically curved surfaces (which is the typical necessary assumption for this theory), in [17], [18] the analysis is made using boundary layer techniques. However, no explicit results can be obtained, as usually with these methods. An early and simplified version of the present result can be found in [16], where incorporating the interface geometry in the asymptotic analysis of a multiscale Darcy-Darcy coupled system is done and a explicit description of the limiting problem is given.
The successful analysis of the present work is due to keeping an interplay between two coordinate systems: the Cartesian and a local one, consistent with the geometry of the interface Γ. While it is convenient to handle the independent variables in Cartesian coordinates, the flow fields in the free fluid region Ω 2 are more manageable when decomposed in normal and tangential directions to the interface (the local system). The a-priori estimates, the properties of weak limits, as well as the structure of the limiting problem will be more easily derived with this double bookkeeping of coordinate systems, rather than trying to leave behind one of them for good. It is therefore a strategic mistake (not a mathematical one, of course) to seek for a transformation flattening out the interface, as it is the usual approach in traces' theory for Sobolev spaces. The proposed method is significantly simpler than other techniques and it is precisely this simplicity which permits to obtain the limiting problem's explicit description for a problem of such complexity, as a multiscale Darcy-Stokes system.

Notation
We shall use standard function spaces (see [20], [1]). For any smooth bounded region G in R N with boundary ∂G , the space of square integrable functions is denoted by L 2 (G ), and the Sobolev space H 1 (G ) consists of those functions in L 2 (G ) for which each of the first-order weak partial derivatives belongs to L 2 (G ). The trace is the continuous linear function γ : H 1 (G ) → L 2 (∂G ) which agrees with restriction to the boundary on smooth functions, i.e., γ(w ) = w ∂G if w ∈ C (cl(G )). Its kernel is The trace space is H 1/2 (∂G ) def = γ(H 1 (G )), the range of γ endowed with the usual norm from the quotient space H 1 (G )/H 1 0 (G ), and we denote by H −1/2 (∂G ) its topological dual. Column vectors and corresponding vector-valued functions will be denoted by boldface symbols, e.g., we denote the product space L 2 (G ) N by L 2 (G ) and the respective N-tuple of Sobolev spaces by H 1 (G ) , furthermore we understand it as a row vector. We shall also use the space H div (G ) of vector functions w ∈ L 2 (G ) whose weak divergence ∇ · w belongs to L 2 (G ). Let n be the unit outward normal vector on ∂G . If w is a vector function on ∂G , we denote its normal component by w n = γ(w) · n, its normal projection by w( n) = w n n. The tangential component is w(tg) = w − w( n). The notation w N , w T indicate respectively, the last component and the first N − 1 components of the vector function w in the canonical basis. For the functions w ∈ H div (G ), there is a normal trace defined on the boundary values, which will be denoted by w · n ∈ H −1/2 (∂G ). For those w ∈ H 1 (G ) this agrees with γ(w) · n. Greek letters are used to denote general second-order tensors. The contraction of two tensors is given by σ : κ = i, j σ ij κ ij . For a tensor-valued function κ on ∂G , we denote the normal component (vector) by κ( n) def = j κ ij n j ∈ R N , and its normal and tangential parts by κ( n) · n = κ( n) n def = i, j κ ij n i n j and κ( n) tg def = κ( n) − κ n n, respectively. For a vector function For a column vector , ... , ∂ ∂x N−1 ; moreover, we regard these operators as row vectors. Finally, ∇ t , ∇ t T denote the same operators written as column vectors, i.e., the operators denoted as column vectors.
Remark 1. It shall be noticed that different notations have been chosen to indicate the first N − 1 components: we use x for a vector as variable x, while we use w T for a vector function w (or the operator ∇ T , ∇ ). This difference in notation will ease keeping track of the involved variables and will not introduce confusion.

Preliminary Results
We close this section recalling some classic results. Lemma 1. Let G ⊂ R N be an open set with Lipschitz boundary, let n be the unit outward normal vector on ∂G . The normal trace operator u ∈ H div (G ) → u · n ∈ H −1/2 (∂G ) is defined by For any g ∈ H −1/2 (∂G ) there exists u ∈ H div (G ) such that u · n = g on ∂G and u H div (G ) ≤ K g H −1/2 (∂G ) , with K depending only on the domain G . In particular, if g belongs to L 2 (∂G ), the function u satisfies the estimate u H div (G ) ≤ K g 0,∂G .
Next we recall a central result to be used in this work Theorem 2. Consider the problem a pair satisfying Here X, Y, X , Y are Hilbert spaces and their corresponding topological duals, F 1 ∈ X , F 2 ∈ Y and the operators A : X → X , B : X → Y , C : Y → Y are linear and continuous. Assume the operators satisfy (i) A is non-negative and X-coercive on ker(B), (iii) C is non-negative and symmetric.
Then, for every F 1 ∈ X and F 2 ∈ Y the problem (2) has a unique solution (x, y) ∈ X × Y, and it satisfies the estimate for a positive constant c depending only on the preceding assumptions on A, B, and C.

Geometric Setting and Formulation of the Problem
In this section we introduce the Darcy-Stokes coupled system analogous to the one presented in [14], for the case when the interface is curved. We begin with the geometric setting

Geometric Setting and Change of Coordinates
We describe here the geometry of the domains to be used in the present work; see Figure 1 (a) for the case N = 2. The disjoint bounded domains Ω 1 and Ω 2 in R N share the common interface, The domain Ω 1 is the porous medium, and Ω 2 is the free fluid region. For simplicity we have assumed that the domain Ω 2 is a cylinder defined by the interface Γ and a small height > 0. It follows that the interface must verify specific requirements for a successful analysis Hypothesis 1. There exists G 0 , G bounded open connected domains G 0 , G ⊂ R N−1 , such that cl(G ) ⊂ G 0 and a C 2 (G 0 ) function ζ : G 0 → R, such that the interface Γ can be described by i.e., Γ is a N − 1 manifold in R N . The domain Ω 2 is described by Remark 2. (i) Observe that the domain G is the orthogonal projection of the open surface Γ ⊆ R N into R N−1 .
(ii) Notice that due to the properties of ζ it must hold that if n = n( x) is the upwards unitary vector, orthogonal to the surface Γ then δ For simplicity of notation in the following we write (iv) For the asymptotic analysis of the coupled system, a domain of reference Ω will have to be fixed, see Figure 1 (b). Therefore we adopt a bijection between domains and account for the changes in the differential operators.

Gradient Operator
Denote by y ∇, x ∇ the gradient operators with respect to the variables y and x respectively. Due to Definition 9 above, a direct computation shows that these operators satisfy the relationship In the block matrix notation above, it is understood that I is the identity matrix in R (N−1)×(N−1) , ∇ T ζ, 0 are vectors in R N−1 and ∂ z = ∂ ∂z . In order to apply these changes to the gradient of a vector function w, we recall the matrix notation reordering we get Here, the operator x D is defined by i.e., x D w ∈ R N×(N−1) and it is introduced to have a more efficient notation. In the next section we address the interface conditions.

Divergence Operator
Observing the diagonal of the matrix in (13) we have Remark 4. The prescript indexes y, x written on the operators above were used only to derive the relation between them , however they will be dropped once the context is clear.

Local vs Global Vector Basis
It shall be seen later on, that the velocities of the channel need to be expressed in terms of an orthonormal basis B, such that the normal vector n belongs to B and the remaining vectors are locally tangent to the interface Γ.
We say the map x → U( x) is a stream line localizer if it is of class C 1 . In the sequel we write it with the following block matrix notation Here, the indexes T and N stands for the first N − 1 components and the last component of the vector field. The indexes tg and n indicate the tangent and normal directions to the interface Γ.
In the following it will be assumed that U is a stream line localizer.
(ii) Notice that by definition U( x) is an orthogonal matrix for all x ∈ G .
Next, we express the velocities in the local basis w 2 in terms of the normal and tangential components, using the following relations Clearly, the relationship between velocities is given by Remark 6. We stress the following observations (i) The procedure above does not modify the dependence of the variables, only the way velocity field are expressed as linear combinations of a convenient (stream line) orthonormal basis.
(ii) The fact that U is a smooth function allows to claim that w 2 tg ∈ H 1 (Ω 2 ) N−1 and w 2 n ∈ H 1 (Ω 2 ). (iii) In order to keep notation as light as possible the dependence of the matrix U as well as the normal and tangential directions n, tg will be omitted whenever is not necessary to show it.
(iv) Notice that given any two flow fields u 2 , w 2 the following isometric identities hold Proposition 3. Let w 2 ∈ H 1 (Ω 2 ) and let w 2 n , w 2 tg be as defined in (19), then 6 (i) Proof. (i) It suffices to observe that the orthogonal matrix U defined in (20) is independent from z.
(ii) Due to (22) we have The last equality holds true since the matrix U( x) is orthogonal at each point x ∈ G , therefore it is an , the result follows.

Interface Conditions and the Strong Form
The interface conditions need to account for stress and mass balance. We start with the stress, to that end we decompose it in the tangential and normal components, the former is handled by the Beavers-Joseph-Saffman (24a) condition and the latter by the classical Robin boundary condition (24b), this gives In the expression (24a) above, 2 is a scaling factor destined to balance out the geometric singularity introduced by the thinness of the channel. In addition, the coefficient α ≥ 0 in (24b) is the fluid entry resistance.
Next, recall that the stress satisfies σ 2 = 2 µ E(v 2 ) (where the scale is introduced according to the thinness of the channel) and that ∇ · v 2 = 0 (since the system is conservative); then we have Replacing in the equations (24) we have the following set of interface conditions The condition (25c) states the fluid flow (or mass) balance.
With the previous considerations, the Darcy-Stokes coupled system in terms of velocity v and pressure p is given by Here, equations (26a), (26b) correspond to the Darcy flow filtration through the porous medium, while equations (26c) and (26d) stand for the Stokes free flow. Finally, we adopt the following boundary conditions The system of equations (26), (27) and (25) constitute the strong form of the Darcy-Stokes coupled system.

Remark 7.
(i) For a detailed exposition on the system's adopted scaling namely, the fluid stress tensor σ 2 = 2 µ E(v 2 ) and the Beavers-Joseph-Saffman condition (24a), together with the formal asymptotic analysis we refer to [15].

Weak Variational Formulation and a Reference Domain
In this section we present the weak variational formulation of the problem defined by the system of equations (26), (27) and (25), on the domain Ω , next, we rescale Ω 2 to get a uniform domain of reference. We begin defining the function spaces where the problem is modeled Let Ω, Ω 1 , Ω 2 , Γ be as introduced in Section 2.1; in particular Ω 2 and Γ satisfy Hypothesis 1. Define the spaces endowed with their respective natural norms. Moreover, for = 1 we simply write X, X 2 and Y.
In order to attain well-posedness of the problem, the following hypothesis is adopted.
Hypothesis 2. It will be assumed that µ > 0 and the coefficients β and α are nonnegative and bounded almost everywhere. Moreover, the tensor Q is elliptic, i.e., there exists a C Q > 0 such that Proof. (i) See Proposition 3 in [14], we simply highlight that the term Ω2 2 β √ Q v 2 (tg) · w 2 (tg) dS has been replaced by Ω2 2 β √ Q v 2 tg · w 2 tg dS, due to the isometric identities (21).
(ii) See Theorem 6 in [14]. The technique identifies the operators A, B, C in the variational statements (29a) and (29b), then it verifies that these operators satisfy the hypotheses of Theorem 2; the result delivers well-posedness.
(iii) A direct substitution of the expressions (14) and (16) in the statements (29), combined with the definition (15) yields the result. Also notice that the determinant of the matrix in the right hand side of the equation (14) is equal to −1 . Finally, observe that the boundary conditions of space X 2 , defined in (28a) are transformed into the boundary conditions of X 2 because none of them involve derivatives.
Remark 8. In order to prevent heavy notation, from now on, we denote the volume integrals by Ω1 F = Ω1 F dx and Ω2 F = Ω2 F d x dz. We will use the explicit notation Ω2 F d x dz only for those cases where specific calculations are needed. Both notations will be clear from the context.

Asymptotic Analysis
In this section, we present the asymptotic analysis of the problem i.e., we obtain a-priori estimates for the solutions (v , p ) : > 0 , derive weak limits and conclude features about those (velocity and pressure) limits. We start recalling a classical space.

Definition 4.
Let Ω 2 as in Definition 1 and define the Hilbert space endowed with its natural inner product Lemma 5. (i) Let H(∂ z , Ω 2 ) be the space introduced in Definition 4, then the trace map w → w Γ from H(∂ z , Ω 2 ) to L 2 (Γ) is well-defined. Moreover, the following Poincaré-type inequalities hold in this space for all w ∈ H(∂ z , Ω 2 ).
Proof. (i) The proof is a direct application of the fundamental theorem of calculus on the smooth functions C ∞ (Ω 2 ) which is a dense subspace in H(∂ z , Ω 2 ).
(ii) A direct application of equations (32) on each coordinate of w ∈ H(∂ z , Ω 2 ) delivers the result.
(iii) It follows from a direct application of (i) and (ii) on w 2 n , w 2 tg respectively.
Next we show that the sequence of solutions is globally bounded under the following hypotheses.
Hypothesis 3. In the following, it will be assumed that the sequences (f 2, : > 0) ⊆ L 2 (Ω 2 ) and (h 1, : > 0) ⊆ L 2 (Ω 1 ) are bounded, i.e., there exists C > 0 such that Theorem 6 (Global a-priori Estimate). Let [v , p ] ∈ X×Y be the solution to the Problem (30). There exists a constant K > 0 such that Proof. Set w = v in (30a) and ϕ = p in (30b) and add them together. In addition, apply the Cauchy-Bunyakowsky-Schwartz inequality to the right hand side and recall the Hypothesis 2, this gives The mixed terms were canceled out on the diagonal. We focus on the summand involving an integral in the right hand side of the expression above, i.e., The second inequality holds due to Poincaré's inequality, given that p 1, = 0 on ∂Ω 1 − Γ, as stated in Equation (27a). The equality holds due to (26b). The third inequality holds because the tensor Q and the family of sources (h 1, : > 0) ⊂ L 2 (Ω 1 ) are bounded as stated in Hypothesis 2 and Hypothesis 3 (Equation (34)), respectively. Next, we control the L 2 (Ω 2 )-norm of v 2, . Since v 2, ∈ H 1 (Ω 2 ) ⊂ H(∂ z , Ω 2 ), the estimates (33) apply; combining them with (37), the bound (34) (from Hypothesis 3) in Inequality (36) we have Using the equivalence of norms · 1 , · 2 for 5-D vectors yields In the next subsections we use weak convergence arguments to derive the functional setting of the limiting problem, see Figure 2 for the structure of the limiting functions.
Corollary 7 (Convergence of the Velocities). Let [v , p ] ∈ X×Y be the solution to the Problem (30). There exists a subsequence, still denoted (v : > 0) for which the following holds.
(ii) There exist χ ∈ L 2 (Ω 2 ) and v 2 ∈ H 1 (Ω 2 ) such that furthermore, ξ satisfies the interface and boundary conditions Proof. (i) (The proof is identical to part (i) Corollary 11 in [14], we write it here for the sake of completeness.) Due to the global a-priori Estimate (35) there must exist a weakly convergent subsequence v 1 ∈ H div (Ω 1 ) such that (39a) holds only in the weak L 2 (Ω 1 )-sense. Because of the hypothesis 3 and the equation (26c), the sequence (∇ · v 1, : > 0) ⊂ L 2 (Ω 1 ) is bounded. Then, there must exist yet another subsequence, still denoted the same, such that (39a) holds in the weak H div (Ω 1 )-sense. Now recalling that the divergence operator is linear and continuous with respect to the H div -norm the identity (39b) follows.
Also, from the strong convergence in the statement (40a), it follows that v 2 is independent from z i.e., (40c) holds.
Again, from (35) we know that the sequence D v 2, : > 0 is bounded in L 2 (Ω 2 ), recalling the identity (15) we have that the expression is bounded. In the equation above the left hand side and the second summand of the right hand side are bounded in L 2 (Ω 2 ), then we conclude that the first summand of the right hand side is also bounded. Hence, ∇v 2, : > 0 is bounded in L 2 (Ω 2 ) and therefore the sequence v 2, : > 0 is bounded in H 1 (Ω 2 ) consequently, the statement (40b) holds.

Theorem 8 (Convergence of Pressures).
Let v , p ∈ X × Y be the solution of problem (30). There exists a subsequence, still denoted (p : > 0) verifying the following.
Proof. (i) (The proof is identical to part (i) Lemma 12 in [14], we write it here for the sake of completeness.) Due to (26b) and (36) it follows that with C > 0 an adequate positive constant. From (27a), the Poincaré inequality implies there exists a constant C > 0 satisfying Therefore, the sequence (p 1, : > 0) is bounded in H 1 (Ω 1 ) and the convergence statement (44a) follows directly. Again, given that p 1, satisfies the Darcy equation (26b) and that the gradient ∇ is linear and continuous in H 1 (Ω 1 ) the equality Qv 1 + ∇p 1 = 0 in (44b) follows. Finally, since p 1, Ω1−Γ = 0 for every element of the weakly convergent subsequence and the trace map ϕ → ϕ Γ is linear, it follows that p 1 satisfies the boundary condition in (44b).
(51) 14 In the expression above the first summand of the second line needs further analysis, we have Combining (48) with the expression above we conclude Back to the equation (51), the two summands on the left hand side of the first line are bounded by a multiple of φ 0,Ω2 due to (48). The first two summands on the third line are trace terms which are also controlled by a multiple of φ 0,Ω2 due to (48). The third summand on the third line is trivially controlled by φ 0,Ω2 due to its construction. Combining all these observations with (52) we get where C > 0 is a new generic constant. The last inequality holds since all the summands in the parenthesis are bounded due to the estimates (35), (46) and the Hypothesis 3. Taking upper limit when → 0, in the previous expression we get The above holds for any φ ∈ C ∞ 0 (Ω 2 ), then the sequence (p 2, : > 0 ) ⊂ L 2 (Ω 2 ) is bounded and consequently the convergence statement (45) follows.
(iii) From the previous part it is clear that the sequence [p 1, , p 2, ] : > 0 is bounded in L 2 (Ω), therefore p also belongs to L 2 Ω) and the proof is complete.
Remark 9. Notice that the upwards normal vector n orthogonal to the surface Γ is given by the expression and the normal derivative satisfies We use the identities above to identify the dependence of χ, ξ and p 2 , see Figure 2.
Theorem 9. Let χ, ξ be the higher order limiting terms in Corollary 7 (ii) and (iii). Let p 2 be the limit pressure in Ω 2 in Lemma 8 (ii). Then In particular, notice that ∂ z ξ = ∂ z ξ( x).
Proof. Take Φ = 0, ϕ 2 ∈ Y, test (30b) and reordering the summands conveniently, we have Letting ↓ 0 in the expression above we get Since the above holds for all ϕ 2 ∈ L 2 0 (Ω 2 ) it follows that where c is a constant. In the expression above we observe that two out of three terms are independent from z, then it follows that the third summand is also independent from z. Since the vector (−∇ T ζ, 1) is independent from z we conclude that ∂ z ξ = ∂ z ξ( x) this, together with the boundary conditions (41b) yield the second equality in (55a).
In the expression above, the first summand clearly tends to zero when ↓ 0. Therefore, we focus on the second summand All the terms in the left hand side can pass to the limit. Recalling the statement (40a), we conclude Letting ↓ 0 in (57) and considering the equality above we get We develop a simpler expression for the sum of the fourth, fifth and sixth terms Here ∂ n is the normal derivative defined in the identity (54b). We introduce the equality above in (58), this yields Next, we integrate by parts the second summand in the first line, add it to the first summand and recall that ∂ z w 2 = −Ψ by construction, thus In the expression above we develop the surface integrals as integrals over the projection G of Γ on R N−1 , we get Recalling that Ψ( x, z) dz · n on Γ, the equality above transforms in Introducing the latter in (60) we have Since the above holds for all Ψ ∈ C ∞ 0 (Ω 2 ) N , we conclude In order to get the normal balance on the interface we could repeat the previous strategy but using Ψ ∈ C ∞ 0 (Ω 2 ) N such that Ψ = Ψ · n n, i.e. such that it is arallel to the normal direction. This would be equivalent to replace Ψ by Ψ · n n in all the previous equations. Consequently, in order to get the normal balance, it suffices to apply · n |(−∇ T ζ, 1)| 2 to Equation (61); such operation yields: In the expression above the identity (42) has been used. Also notice that all the terms are independent from z, then the equation (55b) follows. Consequently, all the terms but the last in (61) are independent of z, therefore we conclude χ is independent from z. Recalling (42) and (55a) the second equality in (55c) follows and the proof is complete.

The Limiting Problem
In this section we derive the form of the limiting problem and characterize it as a Darcy-Brinkman coupled system, where the Brinkman equation takes place in a N − 1-dimensional manifold of R N . First, we need to introduce some extra hypotheses to complete the analysis.
Definition 5. Let x → U( x) be the matrix map introduced in Definition 2. Define the space X tg ⊆ X 2 by endowed with the H 1 (Ω 2 )-norm.
We have the following result Lemma 10. The space X tg ⊂ X 2 is closed.
Proof. Let w 2 ( ) : ∈ N ⊂ X tg and w 2 ∈ X 2 be such that w 2 − w 2 1,Ω2 → 0. We must show that w 2 ∈ X tg . First notice that the convergence in X 2 implies w 2 − w 2 0,Ω2 → 0. Recalling (20) and the fact that U( x) is orthogonal, we have In the identity above, we notice that w 2 T , w 2 T are convergent in the H 1 -norm and the orthonormal matrix U, has differentiability and boundedness properties. Therefore, we conclude that w 2 tg is convergent in the H 1 -norm, we denote the limit by w 2 tg = w 2 tg x, z . Now take the limit in the expression above in the L 2 -sense; there are no derivatives involved, then we have Observe that the latter expression implicitly states that w 2 tg = w 2 tg ( x). Finally, applying the inverse matrix again we have where the equality is in the L 2 -sense. But we know that w 2 tg ∈ H 1 Ω 2 N−1 , therefore the equality holds in the H 1 -sense too, i.e. X tg is closed as desired.
Next, we use space X tg to determine the limiting problem in the tangential direction.
Lemma 11 (Limiting Tangential Behavior's Variational Statement). Let v 2 be the limit found in Theorem 7 (ii). Then the following weak variational statement is satisfied Proof. Let w 2 ∈ X tg , then 0, w 2 ∈ X, test (30a) and get where the second equality holds by the definition of W n and the last equality holds since p 2 is independent from z (55b). Next, recalling the identities (42), (55a) and (55c), observe that Replacing the last two identities in (73) we conclude that the variational statement (72) holds for every test function in W n . Since the bilinear form of the statement is continuous with respect to the norm · X 0 n , it follows that the statement holds for all element w ∈ X 0 n .

Variational Formulation of the Limit Problem
In this section we give a variational formulation of the limiting problem and prove it is well-posed. We begin characterizing the limit form of the conservation laws Lemma 14 (Mass Conservation in the Limiting Problem). Let v 1 , v 2 be the limits found in Theorem 7 then Proof. Take Φ = ϕ 1 , 0 ∈ Y, test (30b) and let ↓ 0, we have The statement above implies (74a).
Next, we introduce the function spaces of the limiting problem Definition 7. Define the space of velocities by endowed with the natural norm of the space X 0 n X tg . Define the space of pressures by endowed with its natural norm.

Theorem 15 (Limiting Problem Variational Formulation).
Let v 1 , v 2 be the limits found in Theorem 7 and let p 1 , p 2 be the limits found in Theorem 8. Then, they satisfy the following variational problem v, p ∈ X 0 × Y 0 : for all w, Φ ∈ X 0 × Y 0 .
In order to prove that the problem is well-posed we prove continuous dependence of the solution with respect to the data. Test (76a) with v 1 , v 2 and (76b) with p 1 , p 2 , add them together and get Applying the Cauchy-Bunyakowsky-Schwarz inequality to the right hand side of the expression above and recalling that v 2 tg is constant in the z-direction we get Here, the second and third inequality holds because p 1 satisfies respectively the drained boundary conditions (Poincaré's inequality applies) and Darcy's equation as stated in (44a). Finally, the fourth inequality is a new application of the Cauchy-Bunyakowsky-Schwarz inequality for 2-D vectors. Introducing (78) in (77) and recalling Hypothesis (2) on the coefficients Q, α, β and µ we have Recalling (39b), the expression above implies that Next, recalling that w 2 tg is independent from z (see (40c)), it follows that v 2 tg 0,Γ = v 2 tg 0,Ω2 and that ∇v 2 tg 0,Ω2 = ∇ T v 2 tg 0,Ω2 . Therefore (79) yields Next, in order to prove continuous dependence for p 2 recall (61) where it is observed that all the terms are already continuously dependent on the data, then it follows that Finally, in order to prove the uniqueness of the solution, assume there are two solutions, test the problem (76) with its difference and subtract them. We conclude the difference of solutions must satisfy the problem (76) with null forcing terms which implies, due to (80), (81) (82) and (83) that the difference of solutions is equal to zero, i.e. the solution is unique. Since (76) has a solution, which is unique and it continuously depend on the data, it follows that the problem is well-posed. Proof. It suffices to observe that due to Hypothesis 4 the limiting problem (76) has unique forcing terms. Therefore, any subsequence of the solutions (v , p ) : > 0 would have a weakly convergent subsequence, whose limit is the solution of problem (76) (v, p), which is also unique, due to Theorem 15. Hence, the result follows.

Closing Remarks
We finish the paper highlighting some aspects that were meticulously addressed in [14].

A Mixed Formulation for the Limiting Problem
For an independent well-posedness proof of the problem (76), define the operators and Then, the variational formulation of the problem (76) has the following mixed formulation The proof now follows showing that the hypotheses of Theorem 2 are satisfied; the strategy is completely analogous to that exposed in Lemma 17, Lemma 18 and Theorem 19 in [14].