A FULLY-MIXED FORMULATION IN BANACH SPACES FOR THE COUPLING OF THE STEADY BRINKMAN–FORCHHEIMER AND DOUBLE-DIFFUSION EQUATIONS

. We propose and analyze a new mixed finite element method for the nonlinear problem given by the coupling of the steady Brinkman–Forchheimer and double-diffusion equations. Besides the velocity, temperature, and concentration, our approach introduces the velocity gradient, the pseudostress tensor, and a pair of vectors involving the temperature/concentration, its gradient and the velocity, as further unknowns. As a consequence, we obtain a fully mixed variational formulation presenting a Banach spaces framework in each set of equations. In this way, and differently from the techniques previously developed for this and related coupled problems, no augmentation procedure needs to be incorporated now into the formulation nor into the solvability analysis. The resulting non-augmented scheme is then written equivalently as a fixed-point equation, so that the well-known Banach theorem, combined with classical results on nonlinear monotone operators and Babuˇska–Brezzi’s theory in Banach spaces, are applied to prove the unique solvability of the continuous and discrete systems. Appropriate finite element subspaces satisfying the required discrete inf-sup conditions are specified, and optimal a priori error estimates are derived. Several numerical examples confirm the theoretical rates of convergence and illustrate the performance and flexibility of the method.


Introduction
The phenomenon in which two scalar fields, such as heat and concentration of a solute, affect the density distribution in a fluid-saturated porous medium, referred to as double-diffusive convection, is a challenging multiphysics problem. This model has numerous applications, among which we highlight predicting and controlling processes arising in geophysics, oceanography, chemical engineering, and energy technology, to name a few areas of interest. In particular, some of them include groundwater system in karst aquifers, chemical processing, simulation of bacterial bioconvection and thermohaline circulation problems, convective flow of carbon nanotubes, and propagation of biological fluids (see, e.g. [1,6,9,22,32]). In this regard, we remark that much of the research Keywords and phrases. Brinkman-Forchheimer equations, double-diffusion equations, fully-mixed formulation, fixed point theory, mixed finite element methods, a priori error analysis.
in porous medium has been focused on the use of Darcy's law. However, this constitutive equation becomes unreliable to model the flow of fluids through highly porous media with Reynolds numbers greater than one, as in the above applications. To overcome this limitation, a first alternative is to employ the Brinkman model [8], which describes Stokes flows through a set of obstacles, and therefore can be applied precisely to that kind of media. Another possible option is the Forchheimer law [23], which accounts for faster flows by including a nonlinear inertial term. According to the above, the Brinkman-Forchheimer equation (see, e.g. [12,16]), which combines the advantages of both models, has been used for fast flows in highly porous media. Moreover, this fact has motivated the introduction of the corresponding coupling with the so called double-diffusion equations (a system of advection-diffusion equations), through convective terms and the body force.
To the authors' knowledge, one of the first works analyzing the coupling of the incompressible Brinkman-Forchheimer and double-diffusion equations is [28], where well-posedness and regularity of solution for a velocitypressure-temperature-concentration variational formulation is established by combining the Galerkin method with a smallness data assumption. Later on, the global solvability of the coupling of the unsteady double-diffusive convection system under homogeneous Neumann boundary conditions and a linearized version of the Brinkman-Forchheimer equations, was introduced and analyzed in [29]. In particular, it is proved in [29] that the global solvability in L 2 -spaces holds true for the 3-dimensional case. More recently, in [31] a finite volume method was adopted to solve the coupling of the time-dependent Brinkman-Forchheimer and double-diffusion equations. The focus of this work was on the validity of the Brinkman-Forchheimer model when various combinations of the thermal Rayleigh number, permeability ratio, inclination angle, thermal conductivity and buoyancy ratio are considered. This study allowed the evaluation of the control parameters effect on the flow structure, and heat and mass transfer. Meanwhile, an augmented fully-mixed formulation based on the introduction of the fluid pseudostress tensor, and the pseudoheat and pseudodiffusive vectors (besides the velocity, temperature and concentration fields) was analyzed in [14]. In there, the well-posedness of the problem is achieved by combining a fixed-point strategy, the Lax-Milgram and Banach-Nečas-Babuška theorems, and the well-known Schauder and Banach fixed-point theorems. The corresponding numerical scheme is based on Raviart-Thomas spaces of order ≥ 0 for approximating the pseudostress tensor, as well as the pseudoheat and pseudodiffusive vectors, whereas continuous piecewise polynomials of degree + 1 are employed for the velocity, and piecewise polynomials of degree for the temperature and concentration fields. Optimal a priori error estimates were also derived.
We point out that the augmented formulation introduced in [14], and the consequent use of classical Raviart-Thomas spaces and continuous piecewise polynomials to define the discrete scheme, are originated by the wish of performing the respective solvability analysis of the Brinkman-Forchheimer equations within a Hilbertian framework. However, it is well known that the introduction of additional terms into the formulation, while having some advantages, also leads to much more expensive schemes in terms of complexity and computational implementation. In order to overcome this, in recent years there has arisen an increasing development on Banach spaces-based mixed finite element methods to solve a wide family of single and coupled nonlinear problems in continuum mechanics (see, e.g. [4, 5, 10-12, 15, 18, 19]). This kind of procedures shows two advantages at least: no augmentation is required, and the spaces to which the unknowns belong are the natural ones arising from the application of the Cauchy-Schwarz and Hölder inequalities to the terms resulting from the testing and integration by parts of the equations of the model. As a consequence, simpler and closer to the original physical model formulations are obtained.
According to the above bibliographic discussion, the goal of the present paper is to continue extending the applicability of the aforementioned Banach spaces framework by proposing now a new fully-mixed formulation, without any augmentation procedure, for the coupled problem studied in [14,28]. To this end, we proceed as in [18] and introduce the velocity gradient and pseudostress tensors as auxiliary unknowns, and subsequently eliminate the pressure unknown using the incompressibility condition. In turn, we follow [18,19] and adopt a dualmixed formulation for the double-difussion equations making use of the temperature/concentration gradients and Bernoulli-type vectors as further unknowns. Then, similarly to [17,18], we combine a fixed-point argument, classical results on nonlinear monotone operators, Babuška-Brezzi's theory in Banach spaces, sufficiently small data assumptions, and the well known Banach fixed-point theorem, to establish existence and uniqueness of solution of both the continuous and discrete formulations. In this regard, and since the formulation for the double-diffusion equations is similar to the ones employed in [18,19], our present analysis certainly makes use of the corresponding results available there. In addition, applying an ad-hoc Strang-type lemma in Banach spaces, we are able to derive the corresponding a priori error estimates. Next, employing Raviart-Thomas spaces of order ≥ 0 for approximating the pseudostress tensor and Bernoulli vectors, and discontinuous piecewise polynomials of degree for the velocity, temperature, concentration and its corresponding gradients fields, we prove that the method is convergent with optimal rate. This work is organized as follows. The remainder of this section describes standard notation and functional spaces to be employed throughout the paper. In Section 2, we introduce the model problem and derive the fully-mixed variational formulation in Banach spaces. Next, in Section 3 we establish the well-posedness of this continuous scheme by means of a fixed-point strategy and Banach's fixed-point theorem. The corresponding Galerkin system is introduced and analyzed in Section 4, where the discrete analogue of the theory used in the continuous case is employed to prove existence and uniqueness of solution. In Section 5, an ad-hoc Strang-type lemma in Banach spaces is utilized to derive the corresponding a priori error estimate and the consequent rates of convergence. Finally, the performance of the method is illustrated in Section 6 with several numerical examples in 2D and 3D, which confirm the aforementioned rates.

The model problem
In what follows we consider the model introduced in [28] (see also [14], Sect. 2), which is given by a steady double-diffusive convection system in a fluid saturated porous medium. More precisely, we focus on solving the coupling of the incompressible Brinkman-Forchheimer and double-diffusion equations, which reduces to finding a velocity field u, a pressure field , a temperature field 1 and a concentration field 2 , the latter two defining a vector := ( 1 , 2 ), such that with parameters := D̃︀/ and F := D R 1 , where D stands for the Darcy number,̃︀ the viscosity, the effective viscosity, R 1 the thermal Rayleigh number, R 2 the solute Rayleigh number, and is a real number that can be calculated experimentally. In addition, the external force f is defined by with g representing the potential type gravitational acceleration, 1,r the reference temperature, 2,r the reference concentration of a solute, and is another parameter experimentally valued that can be assumed to be ≥ 1 (see [28], Sect. 2 for details). The spaces to which 1,r and 2,r belong will be specified later on. In turn, the permeability, and the thermal diffusion and concentration diffusion tensors, are denoted by K, Q 1 and Q 2 , respectively, all them belonging to L ∞ (Ω). Moreover, the inverse of K and tensors Q 1 , Q 2 , are uniformly positive definite tensors, which means that there exist positive constants K , Q1 , and Q2 , such that Equations (2.1) are complemented with Dirichlet boundary conditions for the velocity, the temperature, and the concentration fields, that is u = u D , 1 = 1,D , and 2 = 2,D on Γ, (2.4) with given data u D ∈ H 1/2 (Γ), 1,D ∈ H 1/2 (Γ) and 2,D ∈ H 1/2 (Γ). Owing to the incompressibility of the fluid and the Dirichlet boundary condition for u, the datum u D must satisfy the compatibility condition In addition, due to the first equation of (2.1), and in order to guarantee uniqueness of the pressure, this unknown will be sought in the space Next, in order to derive a new fully-mixed formulation for (2.1)-(2.5), and unlike [14], we do not employ any augmentation procedure and simply proceed as in [18] (see also [19]). More precisely, we now introduce as further unknowns the velocity gradient t, the pseudostress tensor , the temperature/concentration gradient︀ t , and suitable auxiliary variables depending oñ︀ t , u, and , all of which are defined, respectively, by t := ∇u, In this way, applying the matrix trace to the tensors t and , and utilizing the incompressibility condition div(u) = 0 in Ω, one arrives at tr(t) = 0 in Ω and Hence, replacing back (2.7) in the second equation of (2.6), we find that the model problem (2.1)-(2.4) can be rewritten, equivalently, as follows: Find (u, t, ) and ( ,̃︀ t , ), ∈ {1, 2}, in suitable spaces to be indicated below, such that where the Dirichlet datum for is certainly given by D := ( 1,D , 2,D ). At this point we stress that, as suggested by (2.7), is eliminated from the present formulation and computed afterwards in terms of by using that identity. This fact justifies the fourth equation in (2.8), which aims to ensure that the resulting does belong to L 2 0 (Ω).

The fully-mixed variational formulation
In this section we follow [18,19] to derive a fully-mixed formulation in a Banach spaces framework for the coupled system given by (2.8). To this end, we first multiply the third equation of (2.8) by a test function v associated with the unknown u, which formally yields (2.9) Then, applying the Hölder and Cauchy-Schwarz inequalities, we find that the Forchheimer term, given by the second expression in (2.9), can be bounded as where ℓ, ∈ (1, +∞) are conjugate to each other, that is 1 ℓ + 1 = 1. In this way, while we could continue our analysis with arbitrary values of ℓ and , and hence with u and v belonging to the Lebesgue spaces L 2ℓ (Ω) and L (Ω), respectively, we prefer for simplicity to make the latter to coincide, that is such that 2ℓ = , which gives ℓ = 3 2 and = 3, so that both u and v belong to L 3 (Ω). Consequently, the fact that L 3 (Ω) is certainly contained in L 2 (Ω) and the uniform boundedness of K guarantee that the first term in (2.9) is bounded as well, whereas for the third and fourth ones to be well-defined we need to impose that div( ) and f ( ) lie in L 3/2 (Ω). Now, given ∈ (1, +∞), we introduce the Banach space which is endowed with the natural norm ‖ ‖ div ;Ω := ‖ ‖ 0,Ω + ‖div( )‖ 0, ;Ω ∀ ∈ H(div ; Ω).

Analysis of the coupled problem
In this section we combine classical results on nonlinear monotone operators and the Babuška-Brezzi theory in Banach spaces, with the Banach fixed-point theorem, to prove the well-posedness of (2.27) under suitable smallness assumptions on the data. To that end we first collect some previous results and notations that will serve for the forthcoming analysis.

Preliminaries
We begin by establishing the following abstract result.
We end this section by collecting the inf-sup conditions for the operators and̃︀ (cf. (2.16) and (2.25)), and by stating some fundamental properties of the operator (w) (cf. (2.24)), whose proofs follow from a slight adaptation of Lemmas 3.3 and 3.4 in [18], respectively, reason why details are omitted.
Lemma 3.2. There exist positive constants and̃︀, such that (Ω) with boundedness constant given by R ‖w‖ 0,3;Ω , and there hold the following additional properties

Well-definedness of the fixed point operator
In this section we show that the uncoupled problems (3.20) and (3.22) are well-posed, which means, equivalently, that S and̃︀ S (cf. (3.19) and (3.21)) are indeed well-defined. We begin with the operator S. To this end, we first observe that, given ∈ L 6 (Ω), the problem (3.20) has the same structure as the one in Theorem 3.1. Therefore, in order to apply this abstract result, we notice that, thanks to the uniform convexity and separability of L (Ω) for ∈ (1, +∞), all the spaces involved in (3.20), that is, L 3 (Ω), L 2 tr (Ω) and H 0 (︀ div 3/2 ; Ω )︀ , share the same properties.
Proof. It follows straightforwardly from the definition of (cf. (2.15)), along with the triangle, Cauchy-Schwarz, and Hölder's inequalities. Further details are omitted.
Now, let us look at the kernel of the operator (cf. (2.16)), that is which, proceeding similarly to equation (3.34) of [18], reduces to The following lemma establishes hypothesis (ii) of Theorem 3.1 for .
Bearing in mind the definition of (cf. (2.15)), and using (2.3), we obtain Hence, thanks to Lemma 2.1, equation (2.1b) of [3] with = 3, there exists 1 (Ω) > 0, depending only on |Ω|, such that ∫︁ which, together with (3.28), yields Next, bounding below the second term on the right hand side of (3.29) by 0, employing the fact that , and using the continuous injection As a corollary of Lemma 3.5, taking in particular We now establish the unique solvability of the nonlinear problem (3.20).
Lemma 3.6. For each ∈ L 6 (Ω), the problem (3.20) has a unique solution (⃗ u, ) := ((u, t), ) ∈ H × H 0 (div 3/2 ; Ω). Moreover, there exists a positive constant S , independent of , such that Proof. Given ∈ L 6 (Ω), we first recall from (3.9), (3.12) and (3.13) that , D and are all linear and bounded. Thus, bearing in mind Lemmas 3.4 and 3.5, and the inf-sup condition of given by (3.15) (cf. Lem. 3.2), a straightforward application of Theorem 3.1, with 1 = 3 and 2 = 2, to problem (3.20) completes the proof. In particular, noting from (2.15) that (0) is the null functional, we get from (3.4) that and hence the a priori estimate (3.2) yields For later use in the paper we note here that, applying (3.3), and using again the bounds (3.12) and (3.13) for ‖ D ‖ and ‖ ‖, respectively, the a priori estimate for the second component of the solution to the problem defining S (cf. (3.20)) reduces to with depending only on ‖i 3 ‖, BF , BF , and . Next, we aim to proving the well-posedness of problem (3.22), or, equivalently, the well-definedness of the operator̃︀ S (cf. (3.21)), for which, following Lemma 3.6 of [18], we first establish the corresponding hypotheses required by the Babuška-Brezzi theory in Banach spaces. In this way, and similarly to (3.26) and equation (3.35) of [18], we first let̃︀ V be the kernel of the operator̃︀ (cf. (2.25)), that is︀ Then thẽ︀ V-ellipticity of the operator̃︀ is stated as follows.
Lemma 3.7. There exists a positive constant̃︀ such that Proof. We proceed as in Lemma 3.2 of [18]. In fact, given ⃗ := ( ,̃︀ r ) ∈̃︀ V, we know from (3.33) that ∇ =̃︀ r and ∈ H 1 0 (Ω). Hence, using the fact that Q is a uniformly positive definite tensor (cf. (2.3)), and resorting to the Poincaré inequality with positive constant , and to the continuous injection 6 of H 1 (Ω) into L 6 (Ω) (see, e.g. [30], Thm. 1.3.4), we obtain which gives (3.34) with̃︀ := We are now in position to provide the announced result. More precisely, denoting we have the following lemma.

Solvability analysis of the fixed-point equation
Having proved the well-posedness of the uncoupled problems (3.20) and (3.22), which ensures that the operators S,̃︀ S and T are well defined, we now aim to establish the existence of a unique fixed-point of the operator T. For this purpose, in what follows we will verify the hypothesis of the Banach fixed-point theorem. We begin by providing suitable conditions under which T maps a ball into itself. Lemma 3.9. Given > 0, let W be the closed ball in L 3 (Ω) with center at the origin and radius , and assume that the data satisfy Proof. Given w ∈ L 3 (Ω), from the definition of T (cf. (3.23)) and the a priori estimate for S (cf. (3.31)), we first obtain Then, using (3.35) to bound in the foregoing inequality, noting that ‖w‖ 0,3;Ω ≤ , and performing some minor algebraic manipulations, we arrive at which, thanks to the assumption (3.41), yields ‖T(w)‖ 0,3;Ω ≤ and ends the proof.
We now aim to prove that the operator T is Lipschitz continuous, for which, according to its definition (cf. (3.23)), it suffices to show that both S and̃︀ S satisfy this property. We begin with the corresponding result for S.
As a consequence of Lemmas 3.10 and 3.11, we provide next the Lipschitz continuity of T. for all w, z ∈ L 3 (Ω).
We are now in position to establish the main result concerning the solvability of (2.27).

The Galerkin scheme
In this section we introduce and analyze the corresponding Galerkin scheme for the fully-mixed formulation (2.27). The solvability of this scheme is addressed following basically the same techniques employed throughout Section 3.

Solvability analysis
We begin by proving that (4.4) is well posed, or equivalently that S ℎ (cf. (4.3)) is well defined. Indeed, we remark in advance that the respective proof, being the discrete analogue of the one of Lemma 3.6, makes use again of the abstract result given by Theorem 3.1. Hence, we first set the discrete kernel of , which is given by Then, following the approach from Section 5 of [18], we now prove the discrete inf-sup condition for and an intermediate result that will be used to show later on the strong monotonicity of on V ℎ .
We now establish the discrete strong monotonicity and continuity properties of (cf. (2.15)).
Proof. According to Lemma 4.2 and the discrete inf-sup condition for provided by (4.10) (cf. Lem. 4.1), the proof follows from a direct application of Theorem 3.1, with 1 = 3 and 2 = 2, to the discrete setting represented by (4.4). In particular, the a priori bound (4.16) is consequence of the abstract estimate (3.2) applied to (4.4), which makes use of the bounds for D and ℎ given by (3.12) and (3.13).
We remark here that, proceeding similarly to the derivation of (3.32), we obtain with ,d depending only on BF , BF,d , and d . Next, we aim to show that the discrete operator̃︀ S ℎ is well defined. To this end, we now let̃︀ V ℎ be the discrete kernel of̃︀, that is̃︀ Thus, we can establish a preliminary lemma, whose proof follows almost verbatim the one of Lemma 4.1 (see also [5], Lem. 4.2).
Lemma 4.4. There exist positive constants̃︀ d and̃︀ d such that

18)
and The discrete analogue of Lemma 3.8 is established next.
Lemma 4.5. For each w ℎ ∈ H u ℎ , and ∈ {1, 2}, problem (4.6) has a unique solution ( ⃗ ,ℎ , ,ℎ ) = (( ,ℎ ,̃︀ t ,ℎ ), ,ℎ ) ∈̃︀ H ℎ × H ℎ . Moreover, there exists a positive constant̃︀ S,d , independent of w ℎ , such that (4.20) Proof. We proceed as in Lemma 3.8. In fact, given w ℎ ∈ H u ℎ , we first recall from (3.36) and (3.38) that (w ℎ ) is bounded. Then, given ⃗ ,ℎ := ( ,ℎ ,̃︀ r ,ℎ ) ∈̃︀ V ℎ , we easily deduce from (2.3), (4.19), and simple algebraic manipulations, that which, together with the fact that On the other hand, we notice that, following the same arguments yielding (3.40), we are able to show that In what follows we analyze the fixed-point equation (4.8). We begin with the discrete version of Lemma 3.9, whose proof, being a simple translation of the arguments proving that lemma, is omitted. Lemma 4.6. Given > 0, let W ℎ be the closed ball in H u ℎ with center at the origin and radius , and assume that the data satisfy Next, we address the discrete counterparts of Lemmas 3.10 and 3.11, whose proofs, being almost verbatim of the continuous ones, are omitted. We just remark that Lemma 4.7 below is derived using the strong monotonicity of on V ℎ (cf. (4.14)), whereas thẽ︀ V ℎ -ellipticity of̃︀ (cf. (4.21)) and properties (3.17) and (3.18) are employed to obtain Lemma 4.8. Thus, we simply state the corresponding results as follows.
Lemma 4.7. Let BF,d be given by (4.14). Then, there holds for all w ℎ , z ℎ ∈ H u ℎ .
As a straightforward consequence of Lemmas 4.7 and 4.8, we now state the Lipschitz-continuity of the operator T ℎ (cf. Lem. 3.12). Lemma 4.9. Let us define T,d := −1 BF,d̃︀ S,d , with BF,d and̃︀ S,d satisfying (4.14) and (4.24), respectively. Then, there holds for all w ℎ , z ℎ ∈ H u ℎ . We are now in position of establishing the well-posedness of (4.2).
Theorem 4.10. Given > 0, let W ℎ be the closed ball in H u ℎ with center at the origin and radius , and assume that the data satisfy (4.23) and Then the operator T ℎ has a unique fixed point u ℎ ∈ W ℎ . Equivalently, the coupled problem (4.2) has a unique  Proof. It follows similarly to the proof of Theorem 3.13. Indeed, we first notice from Lemma 4.6 that T ℎ maps the ball W ℎ into itself. Next, it is easy to see from (4.25) (cf. Lem. 4.9) and (4.26) that T ℎ is a contraction, and hence the existence and uniqueness results follow from the Banach fixed-point theorem. In addition, it is clear that the estimates (4.29) and (4.30) follow straightforwardly from (4.20) and (4.22), respectively, whereas combining (4.16) (respectively (4.17)) with (4.20) we obtain (4.27) (respectively (4.28)), which ends the proof.
In order to apply Lemma 5.1, we now observe that the problems (2.27) and (4.2) can be rewritten as two pairs of corresponding continuous and discrete formulations of the type defined by (5.1) and (5.2), namely where the operators (u) and (u ℎ ) are defined as in (3.36). The following lemma provides a preliminary estimate for the error ‖(⃗ u, ) − (⃗ u ℎ , ℎ )‖.
Lemma 5.2. There exists a positive constant̂︀ ST ( ), independent of ℎ, such that Proof. We begin by observing that the continuous and discrete systems of (5.4) satisfy the hypotheses of In turn, proceeding as in (3.46), we get Finally, replacing (5.8) back into (5.7), and using the fact that u ∈ W and u ℎ ∈ W ℎ , we readily obtain (5.6) witĥ︀ ST ( ) := ST (1 + 2 )(1 + ), which ends the proof.
Theorem 5.4. Given > 0, assume that the datum D satisfy Then, there exists a positive constant , independent of ℎ, but depending on , BF , BF,d , d , R , ,d ,̃︀ d , ‖Q ‖ 0,∞;Ω , ‖g‖ 0,Ω , ∈ {1, 2}, and the datum D , such that At this point we remark that (3.49), (4.26), and (5.12) share a similar structure holding for a given > 0 and sufficiently small datum D . However, these assumptions are not fully comparable since T , T,d , and ST ( ), being defined in terms of the unknown constants ‖ 6 ‖, ‖i 3 ‖, , d , d , and̃︀ d , are not explicitly computable. As a consequence, we are not able to check in practice whether the examples to be considered below in Section 6 satisfy those hypotheses. Nevertheless, the numerical results reported there confirm the good performance of the method as well as the predicted rates of convergence, which suggests that the aforementioned constraints on the data are more technical issues of the analysis rather than limitations of the applicability of the numerical scheme. In addition, we stress that they do not necessarily have a physical meaning, but only constitute sufficient conditions guaranteeing that problems (2.27) and (4.2) are well-posed, and that the a priori error estimate derived in Theorem 5.4 holds, respectively. In turn, it is important to highlight that (3.49), (4.26), and (5.12) are less restrictive than their counterparts for the augmented fully-mixed formulation proposed in [14], since they only require assumptions on D instead of on u D , D , and r , as in equations (3.46), (4.22), and Theorem 5.4 of [14].

Numerical results
In this section we present three examples illustrating the performance of the fully-mixed finite element method (4.2) on a set of quasi-uniform triangulations of the respective domains, and considering the finite element subspaces defined by (4.1) (cf. Sect. 4.1). In what follows, we refer to the corresponding sets of finite element subspaces generated by = 0 and = 1, as simply P 0 − P 0 − RT 0 − P 0 − P 0 − RT 0 and P 1 − P 1 − RT 1 − P 1 − P 1 − RT 1 , respectively. Our implementation is based on a FreeFem + + code [27], in conjunction with the direct linear solver UMFPACK [20]. A Newton-Raphson algorithm with a fixed tolerance tol = 1E − 6 is used for the resolution of the nonlinear problem (4.2). As usual, the iterative method is finished when the relative error between two consecutive iterations of the complete coefficient vector, namely coeff +1 and coeff , is sufficiently small, that is, where ‖·‖ stands for the usual Euclidean norm in R DOF with DOF denoting the total number of degrees of freedom defining the finite element subspaces H u ℎ , H t ℎ , H ℎ , H ℎ , H̃︀ t ℎ , and H ℎ (cf. (4.1)). We now introduce some additional notation. The individual errors are denoted by: where ℎ stands for the post-processed pressure suggested by the identity (2.7), that is It follows that which shows that the rate of convergence for is at least the one for , which is indeed confirmed below by the numerical results reported. Next, as usual, for each ⋆ ∈ {︁ u, t, , , ,̃︀ t , }︁ we let r(⋆) be the experimental rate of convergence given by where ℎ and̂︀ ℎ denote two consecutive meshsizes with errors e and̂︀ e, respectively. The examples to be considered in this section are described next. Similarly to Section 6 of [14], in all them we take for sake of simplicity = 1, = 1, R 1 = 1, R 2 = 1 and r = (0, 0). In turn, in the first two examples the tensors K, Q 1 , and Q 2 are taken as the identity matrix I, which satisfy (2.3). In addition, the mean value of tr( ℎ ) over Ω is fixed via a Lagrange multiplier strategy (adding one row and one column to the matrix system that solves (4.4) for u ℎ , t ℎ , and ℎ ).

Example 1: 2D domain with different values of the parameter F
In this example we replicate the one from Section 6, Example 1 of [14]. More precisely, we corroborate the rates of convergence in a two-dimensional domain and also study the performance of the numerical method with respect to the number of Newton iterations required to achieve certain tolerance when different values of the parameter F are given. The domain is the square Ω = (−1, 1) 2 . We consider the potential type gravitational acceleration g = (0, −1) t , and choose the data f (cf. (2.2)) such that the exact solution is given by Table 1. Example 1: Number of degrees of freedom, meshsizes, Newton iteration count, errors, and rates of convergence for the fully-mixed P 0 − P 0 − RT 0 − P 0 − P 0 − RT 0 approximation for the coupling of the Brinkman-Forchheimer and double-diffusion equations with F = 10.
The model problem is then complemented with the appropriate Dirichlet boundary conditions. Tables 1 and 2 show the convergence history for a sequence of quasi-uniform mesh refinements, including the number of Newton iterations when F = 10. Notice that we are able not only to approximate the original unknowns but also the pressure field through the formula (6.1). The results confirm that the optimal rates of convergence (ℎ +1 ) predicted by Theorem 5.5 are attained for = 0, 1. The Newton method exhibits a behavior independent of the meshsize, converging in five iterations in all cases. In Figure 1 we display the solution obtained with the fully-mixed P 1 − P 1 − RT 1 − P 1 − P 1 − RT 1 approximation with meshsize ℎ = 0.0284 and 39 102 triangle elements (actually representing 2 074 454 DOF). On the other hand, in Table 3 we show the behaviour of the iterative method as a function of the parameter F ∈ {10 0 , 10 1 , 10 2 , 10 3 , 10 4 , 10 5 }, considering polynomial degree = 0, different meshsizes ℎ, and a tolerance tol = 1E − 06. In this way, here we illustrate that the inertial term F |u|u is well handled by the mixed finite element method (4.2), and that the latter evidences a robust behavior with respect to the parameter F. In fact, only 9 Newton iterations are required to converge in the most challenging case, namely F = 10 5 .

Example 2: Convergence against smooth exact solutions in a 3D domain
We now replicate ( [14], Sect. 6, Example 2). More precisely, we consider the cube domain Ω = (0, 1) 3 and the exact solution: Similarly to the first example, we consider F = 10 and g = (0, 0, −1) t , whereas the data f is computed from (2.2) using the above solution. The numerical solutions are shown in Figure 2, which were built using the fully-mixed approximation with meshsize ℎ = 0.0643 and 63 888 tetrahedral elements (actually representing 1 867 272 DOF). The convergence history for a set of quasi-uniform mesh refinements using = 0 is shown in Table 4. Again, the mixed finite element method converges optimally with order (ℎ), as it was proved by Theorem 5.5.

Example 3: Flow through porous media with channel network
This last example is inspired by Section 5.2.4 of [2], which, similarly to Section 6, Example 3 of [14], focuses on flow through porous media with channel network. To this end, we consider the square domain Ω = (−1, 1) 2 with an internal channel network denoted as Ω c (see the first plot of Fig. 3 below), and boundary Γ, whose left, right, upper and lower parts are given by Γ left = {−1} × (−1, 1), Γ right = {1} × (−1, 1), Γ top = (−1, 1) × {1}, and Γ bottom = (−1, 1) × {−1}, respectively. We consider the coupling of the Brinkman-Forchheimer and doublediffusion equations (2.8) in the whole domain Ω with Q 1 = 0.5 I and Q 2 = 0.125 I, but with different values of the parameters F and K = I for the interior and the exterior of the channel, that is, The parameter choice corresponds to increased inertial effect (F = 10) in the channel and a high permeability ( = 1), compared to reduced inertial effect (F = 1) in the porous medium and low permeability ( = 0.001).
In particular, the first row of boundary equations corresponds to inflow on the left boundary and zero stress outflow on the rest of the boundary. We point out that, differently from Section 6, Example 3 of [14], Dirichlet boundary conditions for temperature and concentration are assumed on the top and bottom of the domain instead of on the left and right sides of Ω as in [14]. We also note that, using similar arguments to those  employed in [15], we are able to extended our analysis to the present case of mixed boundary conditions for the double-diffusion equations. In Figure 3, we display the computed magnitude of the velocity, velocity gradient, pseudostress tensor, and gradients of the temperature and concentration, and the temperature and concentration fields, which were built using the fully-mixed P 0 − P 0 − RT 0 − P 0 − P 0 − RT 0 approximation on a mesh with 27 287 triangle elements (actually representing 475 313 DOF). Similarly to [14], faster flow through the channel network, with a significant velocity gradient across the interface between the porous medium and the channel, are observed here. In addition, the magnitude of the pseudostress tensor is more diffused, since it includes the pressure field. In turn, the temperature and concentration are zero on the top of the domain and go increasing towards the bottom of it, which is consistent with the behavior observed in the magnitude of the temperature and concentration gradients. According to the above, we stress that the fully-mixed approach that    ; computed magnitude of the velocity gradient and pseudostress tensor, and temperature field (middle plots); magnitude of the temperature gradient, concentration field, and magnitude of the concentration gradient (bottom plots).
we have proposed for the coupling of the Brinkman-Forchheimer and double-diffusion equations has the ability to handle heterogeneous media using spatially varying parameters. Moreover, while this example is certainly the most challenging one, due to the strong jump discontinuity of the parameters across the two regions, we highlight that the numerical method (4.2) was able to handle it very efficiently. We notice that the mesh used in this example was built by considering an appropriate refinement around the interface that couples the porous medium with the channel network. Nevertheless, this refinement can be automatized by employing a suitable a posteriori error indicator that captures the aforementioned discontinuity of the parameters. The corresponding a posteriori error analysis and numerical implementation will be addressed in a future work.
We end this paper with a comparison between the present approach and the one from [14] in terms of the corresponding DOF involved and the number of unknowns that they approximate. Indeed, while at first glance the fully-mixed method (4.2) seems more expensive than its augmented mixed counterpart from [14], we stress that the increase observed in the number of unknowns of the Galerkin scheme (4.2) with respect to that from [14] for the same mesh, is due to the fact that, differently from the latter, the former provides direct approximations to three additional variables of physical interest as well, namely the velocity gradient tensor t, the temperature gradient vector̃︀ t 1 , and the concentration gradient vector̃︀ t 2 , in addition to yielding the possibility of employing a post-processing formula to recover the pressure (6.1). Moreover, these four further approximations hold with the same rate of convergence of the remaining variables. However, if one wanted to use the method from [14] to approximate the aforementioned extra unknowns, then one would need to employ numerical differentiation, which, as we know, leads to loss of accuracy of the respective computations.