Robust preconditioning for coupled Stokes-Darcy problems with the Darcy problem in primal form

The coupled Darcy-Stokes problem is widely used for modeling fluid transport in physical systems consisting of a porous part and a free part. In this work we consider preconditioners for monolitic solution algorithms of the coupled Darcy-Stokes problem, where the Darcy problem is in primal form. We employ the operator preconditioning framework and utilize a fractional solver at the interface between the problems to obtain order optimal schemes that are robust with respect to the material parameters, i.e. the permeability, viscosity and Beavers-Joseph-Saffman parameter. Our approach is similar to our earlier work, but since the Darcy problem is in primal form, the mass conservation at the interface introduces some challenges. These challenges will be specifically addressed in this paper. Numerical experiments illustrating the performance are provided. The preconditioner is posed in non-standard Sobolev spaces which may be perceived as an obstacle for its use in applications. However, we detail the implementational aspects and show that the preconditioner is quite feasible to realize in practice.


Introduction
Let Ω = Ω f ∪Ω p , where Ω f is the domain of the viscous flow, Ω p is the domain of the porous media and Γ their common interface. Further let the domain boundaries be decomposed as ∂Ω f = Γ ∪ ∂Ω f,D ∪ ∂Ω f,N and ∂Ω p = Γ ∪ ∂Ω p,D ∪ ∂Ω p,N , where subscripts D, N signify respectively that Dirichlet and Neumann boundary conditions are prescribed on the part of the boundary. The boundary of Γ, i.e., the intersection of Γ and ∂Ω is denoted by ∂Γ. An illustration is given in Figure 1.
The Stokes problem reads: while the Darcy problem in primal form reads: Here, u f , p f are the unknown velocity and pressure for the Stokes problem (1)- (2) in Ω f , p p is the unknown pressure of the Darcy problem (3) in Ω p . The material parameters are the fluid viscosity µ and the permeability K. Here we shall consider the problem with the Dirichlet boundary conditions f on ∂Ω f,D , p p = p 0 p on ∂Ω p,D and Neumann conditions (µ∇u f − p f I) · n f = h on ∂Ω f,N , ∇p p · n p = h p on ∂Ω p,N , where n f , n p are the outer unit normals of the respective subdomains. In particular we assume that |∂Ω i,D | > 0 and |∂Ω i,N | > 0 for i = p, f . Moreover, the coupled problem must be equipped with interface conditions expressing the continuity of stress as well as mass balance. We postpone their description until we describe the weak formulation of the problem. The discretization of the coupled Darcy-Stokes problem with the Darcy problem in a mixed form is challenging since the Darcy and Stokes problems, respectively, call for different schemes. For example, typical finite element methods for the Darcy problem, like the Raviart-Thomas or Brezzi-Douglas-Marini elements, are not stable for Stokes problem as the discretization of the flux specifically targets the properties of H(div) rather than H 1 which is natural for Stokes discretizations. For this reason, a wide range of methods have been proposed over the last decade that address this particular challenge. For example, new elements robust for both the Darcy and Stokes problem have been proposed in [2,3,4,5,6]. Alternatively, stabilization or modifications of standard methods may be used as in [7,8,9]. In this work we will consider the coupled problem with the Darcy equation in a primal form. Standard elements in both the Darcy and the Stokes domain will be used together with a Lagrange multiplier to couple the unknowns appropriately at the interface.
The well-posedness of the Darcy-Stokes problem coupled together through the use of a Lagrange multiplier is wellknown when the Darcy problem is in mixed form [10,11], where both the continuous setting and various discretizations were proposed. Other solution and discretization algorithms for the coupled problem are presented in e.g. [12,13], see [14,10] for an overview. For the mixed formulation we have, in our previous work [1], developed monolithic solvers that are robust with respect to all material parameters by utilizing fractional solvers on the interface. Here, we continue with the same type of approach, but address the difficulty of the Darcy problem in primal form. We remark that the problem to be studied further is symmetric and includes an explicit variable, the Lagrange multiplier, on Γ. In this respect it differs from the more common primal formulation, which leads to a non-symmetric system to be solved for u f , p f and p p . Well-posedness of the latter problem was established in [14] with efficient solvers proposed and analyzed e.g. in [15,16,17].
An outline of the paper is as follows: Section 2 describes the notation, introduces the symmetric primal Darcy-Stokes problem and illustrates the difficulties in its preconditioning. The main challenge for the solver construction, i.e. the proper posing of the coupling operator, is addressed in Section 3. Parameter robust preconditioners are then established in Section 4.

Preliminaries
Let Ω be a bounded Lipschitz domain in R n , n=2 or 3, and denote its boundary by ∂Ω. We denote by L 2 (Ω) the Lebesgue space of square integrable functions, with the norm u 2 L 2 (Ω) = Ω |u| 2 dx, and by H 1 (Ω) the Sobolev space of functions with first derivative in L 2 (Ω) with norm u 2 H 1 (Ω) = u 2 L 2 (Ω) + ∇u 2 L 2 (Ω) . Note that the spaces are both Hilbert spaces, with the standard inner products. These spaces are defined in the same way when u is a vector field, in which case we will write u in boldface. We also define the subspace H 1 0 (Ω) to be the completion in · H 1 (Ω) of C ∞ 0 (Ω), the space of smooth functions on Ω whose restriction to ∂Ω is zero.
For a Lipschitz domain Ω with Γ ⊂ ∂Ω, we can define a trace operator T by T u = u| Γ for smooth u. This can be extended to a bounded, surjective and right-invertible operator H 1 (Ω) → H 1 2 (Γ) (cf. e.g. [18]), where the space H 1 2 (Γ) will be defined later. Given a subset ∂Ω D of ∂Ω, we let H 1 0,∂Ω D (Ω), or for readability just H 1 0,D (Ω), be the subspace of H 1 (Ω) for which the restriction to ∂Ω D is zero, where the restriction is defined in terms of the trace operator. Typically, ∂Ω D will be the subset of ∂Ω on which Dirichlet conditions are prescribed. We also define the semi-norm L 2 τ (Γ) on H 1 (Ω) to be the L 2 (Γ) norm of the tangential component of u at Γ. In 2D, this is just u| Γ · τ L 2 (Γ) where τ is a tangent unit vector, while in 3D it is more conveniently written as u| Γ − (u| Γ · n)n L 2 (Γ) .
For any inner product space X, we let (·, ·) X denote its inner product. When X = L 2 (Ω), we will omit the subscript if there is no cause for confusion. We write the space of continuous linear operators from X to Y as L(X, Y ), or just as L(X) if Y = X. For any two Sobolev spaces X, Y both contained in a common ambient space, we define the intersection and sum spaces X ∩ Y and X + Y in terms of the norms For any c > 0, we define the scaled space cX to be just X as a set, but with the inner product (u, v) X = c(u, v) X . Its norm is trivially equivalent to · X , but because the equivalence constant depends on c, the distinction between the two norms becomes important when we need to establish the independence of bounds with respect to problem parameters.
We define the fractional space H s (Γ) following [19]. Let S ∈ L(H 1 (Γ)) be the operator such that (Su, v) H 1 = (S(I − ∆)u, v) = (u, v) L 2 for all v ∈ H 1 (Γ). We can then find a basis of H 1 (Γ) of orthonormal eigenfunctions e i of S with eigenvalues λ i > 0. Writing u = i c i e i in this basis, we define the norm u 2 H s (Γ) = c 2 i λ −s i for any s ∈ [−1, 1]. Further, let the space H s (Γ) be the completion of C ∞ (Γ) with respect to · H s (Γ) . We also define the space H s 00 (Γ) in the same manner, except that we then apply Dirichlet boundary conditions by choosing S in L H 1 0 (Γ) . Furthermore, H s 00 (Γ) is the completion of C ∞ 0 (Γ) rather than C ∞ (Γ). For the sake of completeness we review here the construction of a matrix realization of fractional operators given in [19]. To this end let V h ⊂ H 1 (Γ), n = dim V h be a finite dimensional finite element subspace with basis functions φ i , i = 1, . . . , n and A, M ∈ R n×n be the symmetric positive definite (stiffness and mass) matrices such that In case V h ⊂ H 1 (Γ) and piecewise constant (P0) discretization is used we let where N is a set of all the facets of the finite element mesh. Further the (facet) average and jump operators are defined When u is a vector function, we define the normal trace T n u = u| Γ · n using the trace operator T component-wise. As such T n is a continuous map H 1 (Ω) → H 1 2 (Γ). Moreover, we let T t be the tangential trace operator. We remark that in 2D and 3D the operator maps to scalar, respectively vector fields. The normal derivative, ∂ n u = ∇u · n| Γ , is more challenging to define properly in this context. Let us therefore briefly sketch an approach, which at least in the authors' opinion at first glance seems like a natural starting point. However, as we will show, the approach does not yield robust preconditioners in our context. First, notice that if we impose additional regularity on u and require that ∆u ∈ L 2 then ∂ n is well defined. In detail, let w ∈ H 1/2 (∂Ω) and E : H 1/2 (∂Ω) → H 1 (Ω) be a (harmonic) extension operator. Then ∂ n u clearly lies in H −1/2 (∂Ω) because This extra regularity assumption is, however, hard to express in the operator preconditioning framework. In particular, to the author's knowledge, there are no standard finite elements that would enable us to exploit the extra regularity. A possible approach could be NURBS [20] or C 1 discretizations developed for fourth order problems. However, the latter often show poor performance for second order problems [21].
Alternatively, we may attempt to define ∂ n as a composition of the first order derivative operator, ∇, with the 1/2 order normal trace operator T n . The composition ∂ n could then be expected to be a 3/2 operator ∂ n : From an operator preconditioning point of view, this would be feasible to realize, as we will see below. However, as we will demonstrate, robustness will not be obtained if we realize ∂ n as a 3/2 operator. In fact, robustness is only obtained if ∂ n is a first order operator, ∂ n : H 1 (Ω) → L 2 (∂Ω). We remark here that while the operator in a continuous setting is ∂ n : H 1 (Ω) → L 2 (∂Ω), in the discrete setting we will include a scaling parameter, i.e. the mesh size, because we use the finite element method. To see that this is reasonable, notice that for finite elements, the mass matrix, as representation of the identity, is differently scaled in different dimensions. In Example 3.2 we detail the scaling in a simplified example.
In order to demonstrate why posing the ∂ n operator properly is required, let us now formulate the coupled Darcy-Stokes problem, where the Darcy problem is in primal form. As a starting point, let the Lagrangian of the coupled problem be, Note that the sign of p f has been changed from (1). Here, the Lagrange multiplier λ in Γ (T n u f − K∂ n p p )λ ds is used to ensure mass conservation, while the extra term Γ 1 2 D(u f · τ ) 2 ds, where D = α BJS µ K , corresponds to the Beavers-Joseph-Saffman condition [22].
The corresponding weak formulation is obtained by the first order optimality conditions of the Lagrangian, that is; where the bilinear forms a, b are defined as We shall refer to (4) as the (primal) Darcy-Stokes problem. Note that the resulting formulation is symmetric. While appropriate function spaces are readily available for u f , p p , p f and their corresponding test functions, it is less clear what the appropriate requirements are for w and λ. This will be addressed below.
Example 2.1. Preconditioner for coupled Darcy-Stokes problem assuming ∂ n : H 1 → H −1/2 . Let us assume that ∂ n is a 3/2 operator so that K∂ n p p ∈ 1 In turn, we consider the following weak formulation: Find The coefficient matrix associated with (4) reads Assuming that the proposed spaces indeed lead to well-posed operator A, the operator preconditioning framework [23] yields as a preconditioner the Riesz mapping In order to test the preconditioner, we solve problem (6) on Figure 1. The mesh is a uniform triangular mesh, consisting of 4N 2 equally sized isosceles triangles. To discretize (4), we use lowest order (P2-P1) Taylor-Hood elements for the Stokes velocity and pressure, while piecewise quadratic elements (P2) were used for the Darcy pressure and piecewise constant elements (P0) for the Lagrange multiplier. Discretization is carried out in the FEniCS library [24], with coupling maps between the interface and domains and the fractional Laplacians being implemented by the extension FEniCS ii [25].
Approximation of the preconditioner (8) is constructed by using single sweep of V -cycle of algebraic multigrid Boomer-AMG from the Hypre library [26] for all the blocks except for the interface block, which is inverted exactly. Starting from a random initial vector, we count the number of iterations required to solve the preconditioned linear system using the MINRES solver from the PETSc library [27] with convergence criterion based on relative tolerance of 10 −8 and absolute tolerance of 10 −10 . Additionally, the condition numbers of B −1 A are computed using an iterative solver from the SLEPc library [28]. In the condition number computations the operator B is computed exactly, that is, all the blocks are inverted by LU. We remark that the solver setup should be used also in the subsequent examples.
The results of the experiment are plotted in Figure 2. By the failure of the iteration counts to stabilize, we see that using 1 µ (I + ∆) −1/2 + K (I + ∆) 1/2 as multiplier space does not lead to a robust preconditioner over the whole parameter range. Note, however, that in the regime where µ is significantly smaller than K (i.e. the lower left region of the plots in Figure 2), iteration counts and condition numbers appear to be stable as the mesh is refined. In this regime, the norm of the multiplier space is dominated by the part from 1 √ µ H −1/2 , which is determined by posing of the trace operator. This suggests that the choice of √ KH 1/2 , i.e. wrong posing of the ∂ n operator, is responsible for the lack of boundedness.

Approximating the trace normal gradient operator
A crucial step in the analysis of the Darcy-Stokes problem will be the mapping properties of the operator ∂ n . As a computationally practical choice of space for the Darcy pressure is √ KH 1 , we immediately run into the problem discussed in the preliminaries because ∂ n cannot be defined on all of H 1 . This necessitates either an assumption of extra regularity or an alternative approach.
Motivated by the observation in [29], that in a discrete finite element setting the trace operator is stable as a map L 2 (Ω) → L 2 (∂Ω), we propose an alternative approach to construct the preconditioners. We start off by outlining the construction of an operator ∂ n, : H 1 (Ω p ) → L 2 (Γ) which will be an approximation to ∂ n . Suppose Γ is a sufficiently regular subset of ∂Ω p , and that Γ is of co-dimension 1 in Ω p . The -thick envelope Γ = {y ∈ Ω p , dist(y, Γ) < } is a higher-dimensional approximation of Γ. For any v ∈ H 1 (Ω p ), where φ is a test function in H 1 (Ω p ).
Note that although the integral over Γ is not well-defined for a general v ∈ L 2 (Ω p ), the integral over Γ is. Provided Γ is sufficiently regular and sufficiently small, we assume that there exists a vector field n Γ on Γ which approximates the normal vector n Γ of Γ at Γ. Using n Γ , we further assume that we can define a bounded extension E : L 2 (Γ) → L 2 (Γ ) along n Γ for which Γ w ds ≈ 1 Γ E w dx for any w ∈ L 2 (Γ). Provided n Γ and E can be defined, then for any u ∈ H 1 (Ω p ) we can define ∂ n, u by Γ ∂ n, u · w ds = 1 Γ ∇u · n Γ E w dx for any w ∈ L 2 (Γ), thus defining the required map ∂ n, : H 1 (Ω p ) → L 2 (Γ ) approximating ∂ n . We assume that the resulting operator ∂ n, is both surjective and bounded, with ∂ n, u L 2 (Γ) ≤ C u H 1 (Ωp) , and that ∂ n, has a bounded right inverse. We emphasize that ∂ n, is just an analytical tool constructed for the analysis in the continuous setting and that is not related to the mesh size h. In fact, we can choose far smaller than the mesh size and for any practical purposes in computations we assume that ∂ n, will be practically identical to ∂ n . We summarize the assumption as follows: Assumption 1. Given a sufficiently regular Γ, ∂ n, : H 1 (Ω p ) → L 2 (Γ) is a bounded surjection which approximates ∂ n on the subspace of H 1 on which ∂ n can be defined. Further, ∂ n, has a bounded right inverse.
Although characterizing the conditions under which Assumption 1 holds is beyond the scope of this paper, we motivate the existence of the required constructions E , n in a few simple examples below.
Example 3.1. Let Γ be the y−axis, and Ω p be the positive half-plane. The construction of E , n Γ is then given by n Γ = n Γ = (−1, 0) and for w(y) ∈ C 1 (Γ) we let (E w)(x, y) = w(y). This continuously extends to all of L 2 (Γ). Clearly ∂ n, u → ∂ n u as → 0 for u ∈ C 1 . Given any w(y) ∈ C 0 (Γ), define u by u(x, y) = −xw(y). Then the map w → u continuously extends to a right inverse of ∂ n, , as by linearity ∂ n, u = ∂ n u = w.
Next, suppose Ω p is the unit disk, and Γ its boundary. By parametrizing Γ with e.g. polar coordinates, this case can be effectively translated to the above. n Γ is now the unit radial vector i r , and for any w(θ) ∈ C 1 (Γ), (E w)(r, θ) = w(θ). Again, this definition of E extends to all of L 2 (Γ). Because 1 Γ ∇u · n Γ E w dx = f (r) dr → f (1) as → 0, we again have ∂ n, u → ∂ n u as → 0 for u ∈ C 1 . Analogously to the previous case, a right inverse can be defined by sending any w(θ) ∈ C 0 (Γ) to u(r, θ) = rw(θ).
Before considering the Darcy-Stokes problem, we justify Assumption 1. First we consider a simplified example in order to illustrate how the scaling of mass matrices in different dimensions affect preconditioners constructed via the application of trace operators. Then, in Example 3.3 we construct preconditioners for a Poisson problem with a ∂ n -constraint which is to be enforced by a Lagrange multiplier, cf. the Babuška problem [30] involving the trace operator.
Letting p denote the Lagrange multiplier associated with the boundary constraint, the extrema u ∈ V , p ∈ Q = L 2 (Γ) of the Lagrangian of (10) satisfy the variational problem: Find u ∈ V and p ∈ Q such that The operator of the preconditioned continuous problem then reads where S is to be constructed such that the condition number is bounded in the discretization parameter h. Here we shall consider three constructions. We remark that when using the finite element method, the identity or the mass matrix has eigenvalues such that both the smallest and the largest eigenvalues scale as h d on uniform mesh. First we consider S = I, with eigenvalues ≈ h. Then, following [29], we let S = h −1 I, i.e., a matrix with eigenvalues ≈ 1. Finally, the choice of S = (−∆ + I) −1/2 is included to show that the relevant trace space in (12) is not (by viewing the trace as an order 1/2 operator) H 1/2 so that dual variable would reside in H −1/2 . We remark that the first two operators are in practical computations assembled as weighted mass matrices where the weights for the respective operators are 1 and inverse cell volume. Recalling Preliminaries §2 the matrix representation of the fractional operator is H(1/2).
To compare the three preconditioners, we let Ω be a unit square, Γ = {(x, y) ∈ ∂Ω, x = 0}. Further, the domain shall be discretized uniformly into 4N 2 isosceles triangles with size h = 1/N , see Figure 3. Considering finite element discretization by P2-P1 elements Table 3.1 lists spectral condition numbers of (12). It can be seen that only the S = h −1 I preconditioner leads to results independent of h.
The growth of the condition number in Table 3.1 due to the preconditioner with −1/2 power indeed confirms that H 1/2 is not appropriate in our setting. An attempt to establish the trace space could be based on viewing the trace as an 1/2 operator. Starting from L 2 a formal calculation then leads to the space H −1/2 and H 1/2 as the multipler space. While we do not include here the results for S = (−∆ + I) 1/2 we remark that the condition number behaves practically identically to S = I.   Table 3.2: Condition numbers of (12) with preconditioner using S = h −1 I. Boundedness with different types of triangulations, cf. Figure 3, and discretizations can be observed.
In order to verify that the properties of h −1 I preconditioner are not due to the highly structured mesh, we consider two additional discretizations of Ω shown in Figure 3. In particular, the triangulations are obtained as refinements of the unstructured meshes where in one case the mesh size is uniform while in the other one the mesh is finer close to the multiplier domain Γ. Moreover, using these triangulations, problem (12) shall be discretizated by P2-P1 elements as well P2-P0 elements to provide more evidence for the preconditioner construction. Indeed, Table 3.2 shows that the condition numbers of (12) are bounded irrespective of the underlying mesh and the finite element discretization considered.
With p the Lagrange multiplier associated with ∂ n -constraint (13) leads to a variational problem: Find u ∈ V and p ∈ Q = L 2 (Γ) such that The preconditioned continuous problem then reads Following the preliminaries where ∂ n was regarded as a 3/2 operator we let S = (−∆ + I) 1/2 . Alternatively, S = h −1 I is set following the Assumption 1. Finally S = I is considered. Matrix realization of the S operators shall be identical to Example 3.2. We shall also use the tessellations described in Example 3.2 as well as identical eigenvalue solvers.
To compare the three preconditioners we let Ω be a unit square and Γ = {(x, y) ∈ ∂Ω, x = 0} and we consider first the (Neumann) case where ∂Ω N = {(x, y) ∈ ∂Ω, y = 0 or y = 1}, i.e. where the multiplier domain intersects the part of boundary with Neumann boundary conditions. Using the uniform meshes (marked as (us) Figure in 3) and P2-P1 elements, Table 3.3 shows the spectral condition numbers of (15). As in Example 3.2 only S = h −1 I preconditioner (based on Assumption 1) leads to results independent of h.  Table 3.3: Condition numbers of (15) discretized by P2-P1 elements on uniform refinements of (us) mesh in Figure 3. Boundednes in discretization is obtained only with S = h −1 I.   Table 3.4: Condition numbers of (15) using S = h −1 I preconditioner discretized on uniform refinements of parent meshes in Figure  3 using two element types. Refinement level is indicated by l. (Left) Γ intersects ∂ΩN . (Right) Γ intersects ∂ΩD.
In the context of multiscale problems, compatibility of boundary conditions of the multiplier space and the boundary conditions prescribed on the domain intesecting Γ is known to present an issue, cf. e.g. [11]. Here, we address this problem by considering (15) with |∂ N Ω| = 0, i.e. we let Γ intersect only the Dirichlet boundary. We remark that until this point only intersection with Neumann boundary was considered.
In Table 3.4 the Dirichlet problem is considered with an unmodified h −1 I preconditioner. In particular, with P2-P1 discretization we impose no boundary conditions on the multiplier space. Using this construction the condition numbers can be seen to remain bounded on all the meshes and with both finite element discretizations.
We remark that the h −1 I preconditioner is equally unaffected by the Dirichlet boundary conditions on ∂Ω D = ∂Ω \ Γ in the trace-constrained L 2 projection problem (10) with V = H 1 0,∂Ω D (Ω), cf. Example 3.2, in contrast to the H 1 problems considered in [1], where the appropriate preconditioner was H − 1 2 00 or H − 1 2 depending on whether the interface intersected the Dirichlet boundary or not. We remark that in the continuous setting boundary values have measure zero and this may then be perceived as the L 2 space being the correct one in our discrete setting. Of course, the counterargument in the continuous setting is that then the trace cannot be defined. However, in the discrete setting, this can be done.
Without including the simulation results we comment here that the condition numbers of the Dirichlet problem are practically identical to those presented in Tables 3.1 and 3.2. In addition, with the two preconditioners S = I and S = (−∆ + I) 1/2 on the unstructured meshes a growth of condition numbers with h is observed similar to Table 3.3.
We remark that the stability of the preconditioner h −1 I in Example 3.3 provides numerical evidence for well-posedness of (14), i.e. the Darcy subproblem in the coupled Darcy-Stokes system (4).

Robust Preconditioners for the Darcy-Stokes system
In Example 2.1, we showed that the efficiency of the preconditioner (8) for the primal Darcy-Stokes problem (7) varied substantially with the material parameters even though the Stokes block and the Darcy block were preconditioned with appropriate preconditioners, and argued that the reason was a poor preconditioner at the interface.
In this section we demonstrate that robustness with respect to mesh resolution and variations in material parameters can be obtained by posing the Lagrange multiplier in properly weighted fractional spaces, namely the intersection space No modifications of the velocity or pressure space norms will be required. Our analysis is closely related to [1], and based on Assumption 1 along with an assumption of stability for the Stokes problem. We remark that although Assumption 1 is motivated by the discrete problem, our analysis is carried out in a continuous setting.
We shall prove well-posedness of the coupled Darcy-Stokes problem (4) with spaces We remark that in case Γ intersects only the Dirichlet boundary ∂Ω f,D the space H −1/2 needs to be modified to reflect H 1/2 00 as the appropriate trace space of V f . We refer to [1] for a thorough discussion of the subject. As a prerequisite for the coupled problem to be well-posed, we require that each subproblem is well-posed. For the Stokes subproblem the property has been demonstrated by numerical experiments in [1]. Here we state the result without proof.
are arbitrary. Then we assume that the Stokes problem: Find satisfies the Brezzi conditions and hence has a unique solution (u f , p f , λ) ∈ V S and the following bound holds Here the constant C depends only on Ω f , ∂Ω f,D and Γ.
Corresponding well-posedness of the Darcy problem with ∂ n -constraint was demonstrated numerically for K = 1 in Example 3.3. Here, we analyze the general case.

Lemma 1.
Suppose Ω p , Γ are such that Assumption 1 holds and |∂Ω p,D | > 0. Then for any f ∈ 1 has a unique solution satisfying where C is a constant depending only on Ω p .
. We consider the left-hand side of (3) as an operator where (Ap p , q p ) = K(∇p p , ∇q p ) Ωp and (Bp p , w) = K(∂ n, p p , w) Γ . The statement of the theorem follows from the Brezzi theory [31] once the Brezzi conditions are verified. That is, we must show that A, B are bounded, A is coercive on ker B and that the inf-sup condition inf q∈Q sup v∈V (Bv, q) ≥ β v q holds for some constant β > 0.
Here the boundedness of A and the coercivity on V are evident. For the latter we recall that |∂Ω p,D | > 0 is assumed and invoke the Poincare inequality. Assumption 1 is needed to show the properties of B. Because Hence, the inf-sup condition holds with β = 1 C . By Brezzi theory, Theorem 1 follows, and the problem is well-posed. As in the argument of [1], we note that due to our use of parameter weighted spaces, all constants are in fact independent of the problem parameters, and the operator preconditioner is therefore robust to parameter variations.
Using operator preconditioning and Theorem 1 a suitable preconditioner for the primal Darcy-Stokes problem (4) is a Riesz map with respect to the inner product of W in (16), that is, the operator We remark that all of the components of the preconditioner can be realized in an efficient, order optimal manner with multilevel schemes. In particular, the only non-standard component here is the multilevel scheme for the fractional operator which, however, has been established in [32].
Example 4.1 (Robust Darcy-Stokes preconditioning). We consider the setup from Example 2.1 while using the operator (21) as preconditioner. As before, the leading blocks of the preconditioner are realized using single algebraic multigrid V -cycle. The multiplier block is then assembled using the eigenvalue decomposition and its inverse is computed by a direct solver.
The obtained iteration and condition numbers are plotted in Figure 4. It can be seen that both quantities are bounded in mesh size N as well as the physical parameters µ, κ and α BJS . Remark 1. Below we consider the validity of Assumption 1 in a continuous and discrete setting. Clearly, in a continuous setting it is easy to find a function that violates the assumption. Consider the case where h while Ω and Γ are both unit sized. Further, let u ∈ H 1 (Ω p ) be a function which is zero in Ω\Γ and has a gradient of 1 in Γ . Recalling our definition of the operator ∂ n, : H 1 (Ω p ) → L 2 (Γ ) by Γ ∂ n, u · w ds = 1 Γ ∇u · n Γ E w dx, we see that ∂ n, u ∈ L 2 (Γ ) is the unit constant function whereas u 1 ≈ √ . Hence ∂ n, ≥ ∂n, u L 2 (Γ ) u H 1 (Ωp ) ≈ 1 √ , which is very large for small . Clearly, this function violates Assumption 1.
The above construction of a function that violates the assumption is however clearly not relevant in our discrete setting as these functions are below the resolution of our finite element mesh. Indeed, in our numerical experiments, we use discrete subspaces of H 1 (Ω p ), so that any function whose gradient is nonzero on Γ also has nonzero gradient at distance h from Γ. This means that if is chosen smaller than h, functions like u above which are zero immediately outside of Γ are not admissible.
For a relevant finite element function u h , constructed as above, i.e., such that u h is zero everywhere except having a gradient of 1 on the finite elements with facets on Γ, assuming that h, we have