A new mixed finite element method for the n-dimensional Boussinesq problem with temperature-dependent viscosity

In this paper we propose a new mixed-primal formulation for heat-driven flows with temperature-dependent viscosity modeled by the stationary Boussinesq equations. We analyze the well-posedness of the governing equations in this mathematical structure, for which we employ the Banach fixed-point theorem and the generalized theory of saddle-point problems. The motivation is to overcome a drawback in a recent work by the authors where, in the mixed formulation for the momentum equation, the reciprocal of the viscosity is a pre-factor to a tensor product of velocities; making the analysis quite restrictive, as one needs a given continuous injection that holds only in 2D. We show in this work that by adding both the pseudo-stress and the strain rate tensors as new unknowns in the problem, we get more flexibility in the analysis, covering also the 3D case. The rest of the formulation is based on eliminating the pressure, incorporating augmented Galerkin-type terms in the mixed form of the momentum equation, and defining the normal heat flux as a suitable Lagrange multiplier in a primal formulation for the energy equation. Additionally, the symmetry of the stress is imposed in an ultra-weak sense, and consequently the vorticity tensor is no longer required as part of the unknowns. A finite element method that follows the same setting is then proposed, where we remark that both pressure and vorticity can be recovered from the principal unknowns via postprocessing formulae. The solvability of the discrete problem is analyzed by means of the Brouwer fixed-point theorem, and we derive error estimates in suitable norms. Numerical examples illustrate the performance of the new schem and its use in the simulation of mantle convection, and they also confirm the theoretical rates of convergence.

strain rate tensors as new unknowns in the problem, we get more flexibility in the analysis, covering also the 3D case. The rest of the formulation is based on eliminating the pressure, incorporating augmented Galerkin-type terms in the mixed form of the momentum equation, and defining the normal heat flux as a suitable Lagrange multiplier in a primal formulation for the energy equation. Additionally, the symmetry of the stress is imposed in an ultra-weak sense, and consequently the vorticity tensor is no longer required as part of the unknowns. A finite element method that follows the same setting is then proposed, where we remark that both pressure and vorticity can be recovered from the principal unknowns via postprocessing formulae. The solvability of the discrete problem is analyzed by means of the Brouwer fixed-point theorem, and we derive error estimates in suitable norms. Numerical examples illustrate the performance of the new schem and its use in the simulation of mantle convection, and they also confirm the theoretical rates of convergence.
1. Introduction. The description of a variety of natural phenomena and engineering problems deal with incompressible quasi-Newtonian flows with viscous heating and buoyancy terms (often referred to as natural convection of fluids). Mantle convection with very large viscosities, waves and currents near shorelines, heat transfer in nanoparticle fluids, creeping thermal plumes, stratified oceanic flows, chemical reactors, and many other examples can be invoked in the context of applicative problems, partial differential equations, and in the construction of numerical schemes [17,18,22,26,30,32,33,40]. In particular, finite element methods approximating the solution of these equations have been developed by the mathematical community in the last two decades. For instance, the model with constant coefficients has been addressed via primal approaches in [11,13,19,23], whereas mixed-type schemes have been employed in [20,24,25], and in particular two different formulations based on a dual-mixed approach for the momentum equation, and a primal and mixed-primal one for the energy equation, are proposed in [21].
Here we advocate to the study of well-posedness and stability of a weak formulation (as well as the construction of mixed finite element schemes) for the Boussinesq equations with thermally-dependent viscosity. In this regard, we remark that not many finite element methods that provide analysis for the case of non-constant viscosity are available in the literature (we may basically refer to [2,3,5,34,36,39,[41][42][43], and some of the references therein). For instance, the authors in [36] (see [35] for the continuous analysis) propose an optimally convergent, primal finite element method for the steady-state problem applied to the flow in channels. In turn, [39] (and a stabilized version of their method, which is suggested in [42]) deal with the unsteady problem, where backward Euler discretization is used in time, and conforming finite elements in space. On the other hand, a conforming finite element method for the problem with temperature-dependent viscosity and thermal conductivity is developed in [34]. In this work, Stokes-stable elements for the velocity and pressure, as well as Lagrange elements for the temperature, and discontinuous piecewise polynomials for the normal heat flux through the boundary, are used to define the associated discrete scheme. Similar hypotheses to those in [34] are also required in primal formulations that may be computationally less expensive than a mixed method, but do not approximate directly other variables of physical interest. The user of a primal scheme would usually appeal to numerical differentiation, with the consequent loss of accuracy.
In the Boussinesq approximation it is also assumed that density fluctuations are linearly related to temperature, and so the last term on the momentum balance is the buoyancy contribution that contains a rescaling implying that, at least for the analysis of the problem, the adimensional number Ra/(Pr Re 2 ) is of the order of the gravity magnitude. In addition, and for sake of notational simplicity of the analysis, here we also take a zero reference temperature and assume that other constants can be absorbed by the model adimensionalization.
With respect to the boundary conditions for (1.1), we assume that u D ∈ H 1/2 (Γ) : = [H 1/2 (Γ)] n , ϕ D ∈ H 1/2 (Γ), and that u D verifies the compatibility condition where ν denotes the unit outward normal on Γ. The construction of the mixedprimal formulation in [5] begins with the introduction of the stress and vorticity tensors, respectively, defined as where ⊗ stands for the tensor product operator, I is the identity matrix in R n×n , the symbol t denotes matrix transpose, and Note that if we had Neumann boundary conditions in the fluid, then σν would be prescribed on Γ. Therefore, after eliminating the pressure p, problem (1.1) is rewritten as: Find (σ, u, γ, ϕ) such that where tr denotes the matrix trace, τ d := τ − 1 n tr(τ ) I is the deviatoric tensor of a given τ ∈ L 2 (Ω), and (1.6e) constitutes a uniqueness condition for the pressure. At this point we remark that, because of the division by µ(ϕ) in (1.6a), integration by parts (after multiplication by a test function) is now possible. However, it can be seen that this leads to the usage of a continuous injection from H 1 (Ω) into L 8 (Ω), as required by the following estimate: with C(Ω) > 0 and valid for any ϕ ∈ H 1 (Ω); u, w ∈ H 1 (Ω); τ ∈ L 2 (Ω), and more important, for Ω ⊂ R d , d ∈ {1, 2}, according to the Sobolev embedding theorem (cf., e.g., [37,Theorem 1.3.5]). This estimate is used in several ways throughout [5], at both continuous and discrete levels (see, [5,Lemmas 3.8,4.5 and 5.3]), and its main purpose is to help in the proof of Lipschitz continuity of the fixed-point operator T (respectively its discrete version T h ) that consequently provides well-posedness of the continuous formulation (respectively the Galerkin scheme).
One of the purposes of this work is to derive a new augmented method for (1.1) that remains valid also for the 3D case. To this end, one can look at e.g. [15] where the strain rate tensor e(u) is considered as a new unknown, in addition to either the stress (or pseudostress) and vorticity tensors. This provides more flexibility in the scheme, as it is no longer necessary (nor much convenient) to divide in (1.6a) by the viscosity to set the gradient free. Instead, the decomposition of the velocity gradient into its symmetric and skew-symmetric parts provides an equation to be integrated by parts. However, proceeding a bit differently than in those works, here we do not impose the symmetry of the stress by testing against the whole space of the skew symmetric tensors, but only against a proper subspace of it, thus yielding what we call an ultra-weak imposition of such constraint. As a consequence, we do not require the vorticity as a further unknown (also interpreted as an associated Lagrange multiplier), and hence the resulting mixed scheme has fewer degrees of freedom without affecting the accuracy of the approximate solutions. Moreover, the pressure and the vorticity are easily recovered by postprocessing formulae that provide the same rates of convergence of the main unknowns. In terms of the analysis, we can suitably modify the approach from [5], thus giving rise to an augmented mixedprimal finite element method using discontinuous piecewise polynomial functions of degree ≤ k to approximate the strain rate and normal heat flux, Raviart-Thomas elements of order k for the stress, and Lagrange elements of order k + 1 for the velocity and temperature. We stress that the already described extension to the 3D case, the new way of imposing the symmetry of the stress with the consequent decrease in the number of unknowns, and the availability of the aforementioned postprocessing formulae, constitute the main contributions of the present work. Let us also clarify that rather than proposing a cheap method for any reformulation or simplification of the governing equations, our underlying motivation is to rigorously derive properties of numerical schemes that preserve the mixed structure of the PDE system.
The rest of this work is organized as follows. In Section 2 we rewrite the Boussinesq problem (1.1) considering the strain rate and stress tensors as new variables, to then derive an augmented mixed-primal formulation whose well-posedness is proved by means of a fixed-point approach. Similarly, in Section 3 we provide the corresponding Galerkin scheme and its associated well-posedness, to then, in Section 4, proceed to derive a priori error estimates and state the respective rates of convergence, including those for the postprocessed approximations of pressure and vorticity, when a particular choice of finite element subspaces is made. Finally, to complement our theoretical results, we present in Section 5 a set of numerical examples that serve to confirm the properties of the proposed schemes.
2. The continuous formulation. In order to avoid the division by µ(ϕ) in (1.6a), we have a look at recent work [15], where the authors develop a augmented mixed finite element method for the Navier-Stokes equations with nonlinear viscosity. This approach relies on the definition of the strain rate tensor as a new unknown where L 2 tr (Ω) = s ∈ L 2 (Ω) : s = s t and tr(s) = 0 . In addition, for each v ∈ H 1 (Ω) we let η(v) be the skew-symmetric part of the velocity gradient tensor ∇v, that is which certainly lies in L 2 skew (Ω) (cf. (1.5)). In this way, bearing in mind (2.1), and noting from (1.4) and (2.2) that γ = η(u), the system (1.6) can be rewritten as: Before continuing in the following section with the derivation of the mixedprimal formulation of (2.3), we find it important to remark at this point that it is also possible to account the case of a nonlinear thermal conductivity K in the forthcoming analysis. In particular, within the context of coupled flow-transport and sedimentation-consolidation models, in which the equation for the concentration ϕ reads as the present one for the energy, we may refer to [8, eq. (2.1)] and [9, eq. (2.1)], where diffusion coefficients depending on ∇ϕ and ϕ, respectively, were considered for the corresponding analyses. More recently, a slight modification of the present model (1.1), in which K was allowed to depend on the temperature ϕ, namely in the form κ(ϕ), with κ : R −→ R + , was introduced and analyzed in [4] by using a fully-mixed variational formulation. In this case, the pseudoheat vector p := κ(ϕ)∇ϕ − ϕu and the temperature gradient are introduced as auxiliary unknowns, so that the energy equation is rewritten as ζ = ∇ϕ , p := κ(ϕ)ζ − ϕu , and div(p) = 0 in Ω. Further details on the analysis of such fully-mixed approach are available in [4].

2.
1. An augmented mixed-primal formulation. Multiplying (2.3a) by a test function τ ∈ H(div; Ω), integrating by parts and using the Dirichlet condition (2.3e), we obtain Then, we multiply (2.3b) and (2.3c) by appropriate test functions, imposing at the same time the symmetry of the stress tensor σ in the form We highlight that (2.4) constitutes what we call an ultra-weak imposition of the symmetry of σ since η(v) : v ∈ H 1 (Ω) is a proper subspace of L 2 skew (Ω). This idea has also been applied in [7].
The weak form of the energy conservation equation is recalled next from [5]: where λ := K∇ϕ · ν ∈ H −1/2 (Γ) is a Lagrange multiplier taking care of the nonhomogeneous Dirichlet boundary condition. Notice that, due to the second term in (2.5) and the right hand side of (2.6), the velocity u must live in H 1 (Ω), since appealing to the continuous injection of H 1 (Ω) into L 4 (Ω), there exist positive constants c 1 (Ω) and c 2 (Ω) such that

The fixed-point argument.
A crucial tool in [5] to prove the well-posedness of the continuous and discrete formulations is a technique that decouples the problem into the mixed formulation of the momentum equation and the primal formulation of the energy equation, which further enables us to rewrite the formulation as a fixed-point problem. Hence, we denote H := H 1 (Ω) × H 1 (Ω) and consider in what follows the operator S : H → H defined by where t is the solution of the problem: Find t ∈ H such that In addition, let S : H → H 1 (Ω) be the operator defined by where ϕ is the first component of a solution to: In this way, by introducing the operator T : H → H as we realize that (2.11) can be rewritten as the fixed-point problem: Find (u, ϕ) ∈ H such that T(u, ϕ) = (u, ϕ). (2.23) As in [5], the objective is to use the Banach fixed-point theorem to prove existence and uniqueness of (2.23). We recall that the key difference in the present work with respect to [5] is in the problem that defines the operator S, and therefore, those results associated to the operator S and the primal formulation of the energy equation will be considered here as well, but we only cite them.

Well-posedness of the uncoupled problems.
In what follows, we consider the norms We first recall some results that will be used for ellipticity purposes.
The following results establish sufficient conditions for the operators S and S being well-defined, equivalently, (2.20) and (2.21) being well-posed.
Then, there exists r 0 > 0 such that for each r ∈ (0, r 0 ), the problem (2.20) has a unique solution t := S(w, φ) ∈ H for each (w, φ) ∈ H such that w 1,Ω ≤ r. Moreover, there exists a constant C S > 0, independent of (w, φ) and r, such that there holds Proof. Given (w, φ) ∈ H, we notice from (2.12) that A φ is bilinear. Then, by using the upper bound of the viscosity function, the Cauchy-Schwarz inequality and the trace theorem with constant c 0 (Ω), we see that for any t, s ∈ H, and therefore, there exists a constant C A > 0 depending only on µ 2 , c 0 (Ω) and the stabilization parameters κ j , such that On the other hand, using (2.13) and (2.8), we obtain that for any t := (t, σ, u), s := (s, τ , v) ∈ H there holds which together with (2.25) implies the existence of a positive constant denoted by A + B , independent of (w, φ) such that and then, using the bounds for the viscosity and the Cauchy-Schwarz and Young inequalities, we obtain for any δ 1 , δ 2 > 0 and any s ∈ H that Then, defining the following constants: and using Lemmas 2.1 and 2.2, it is possible to find a constant α(Ω) := min{α 1 , α 4 , α 5 }, independent of (w, φ), such that Combining the foregoing inequality with (2.26), we get that, for any s ∈ H, there holds Therefore, we easily see that thus proving ellipticity for A φ + B w under the requirement (2.29). Finally, the linearity of the functionals F D and F φ is clear, and from (2.16), (2.17), using the Cauchy-Schwarz inequality, the trace estimates in H(div; Ω) and H 1 (Ω), with constants 1 and c 0 (Ω), and the continuous injection from   .24) is satisfied with C S := 2M S α(Ω) , a constant clearly independent of (w, φ). In addition, the fact that α(Ω) and M S are both independent of r guarantees that C S does not depend on r either.
We stress here that, while the viscosity µ could be assumed to be a monotone decreasing function of the temperature φ, this fact by itself would not help to simplify the solvability analysis of (2.20) (which defines the operator S), since φ does not form part of the corresponding unknowns, but it actually constitutes a datum for that problem. Indeed, all what one needs from µ for the proof of Lemma 2.3 are its upper and lower bounds. A similar remark is valid for the Lipschitzcontinuity of S (cf. Lemma 2.6 below) in which the homonimous property of the viscosity plays a key role. As it will become clear later on, the smallness assumptions on the data for the well-posedness of our problem arises from the main properties of the resulting fixed point operator, and particularly from those inherited from S, so that the eventual monotonicity of µ does not help to remove those conditions either.
Lemma 2.4. For each (w, φ) ∈ H there exists a unique pair (ϕ, λ) ∈ H 1 (Ω) × H −1/2 (Γ) solution of the problem (2.21). Moreover, there exists C S > 0, independent of (ϕ, λ) and r, such that Proof. The results comes from a direct application of the Babuška-Brezzi theory (see [5,Lemma 3.6] or more precisely [20,Lemma 3.4]). In particular, the righthand side of (2.31) is obtained after bounding the functionals F w,φ (cf. (2.18)) and G (cf. (2.19)), respectively, whose resulting bounds are independent of r (cf. [20, eq. (3.39)] for F w,φ ). In this way, since the constants yielding the inf-sup condition of b and the ellipticity of a in the kernel of b, are both independent of r, the corresponding continuous dependence result for (2.21) confirms that C S does not depend on r either.
From the previous two lemmas, it is clear that T is well-defined for any pair (w, φ) ∈ W, where which is nothing but the closed ball in H with center (0, 0) and radius r ∈ (0, r 0 ), with r 0 defined in (2.29). We also notice that a particular choice of stabilization parameters is necessary. We begin by selecting the middle points of the ranges for δ 1 , δ 2 , κ 1 and κ 3 , that is, to then pick κ 2 and κ 4 so that α 2 and α 3 can attain the largest value possible, that is We stress here that the foregoing feasible choices of κ j , j ∈ {1, 2, 3, 4}, are explicitly computable since they do not depend on the usually unknown constants c 3 (Ω) and ϑ 0 (Ω) from Lemmas 2.1 and 2.2, respectively, but only on the known constants µ 1 and µ 2 bounding the viscosity (cf. (1.2)).

2.4.
Further-regularity assumption. Although the problem that defines S (i.e. (2.20)) is well-posed, a small further-regularity assumption is needed in order to continue with the analysis. Inspired by [8], we assume that u D ∈ H 1/2+ε (Γ) for some ε ∈ (0, 1) when n = 2, or ε ∈ 1 2 , 1 when n = 3, and that for each with C S (r) > 0 independent of z but depending on the upper bound r of its H 1norm.
2.5. Solvability analysis of the fixed-point equation. Throughout the rest of the paper we assume that the regularity hypotheses of the previous section are valid. We now proceed to directly fulfill the hypotheses of the Banach fixed-point theorem. The following result shows that T can map the ball W into itself.
Lemma 2.5. Consider the closed ball W defined in (2.32) with r ∈ (0, r 0 ) and r 0 as given in (2.29). Suppose the data satisfy Next, we prove some results that will help us to arrive at the Lipschitz continuity of T. Lemma 2.6. Let r ∈ (0, r 0 ) with r 0 as given in (2.29). Then, there exists a positive constant C S , independent of r, such that
Proof. Let (w 1 , φ 1 ), (w 2 , φ 2 ) ∈ H as indicated and let t j := (t j , σ j , u j ) = S(w j , φ j ) ∈ H, j ∈ {1, 2} be the corresponding solutions of (2.20). Then, adding and subtracting the equality we find that Thus, using the ellipticity of A φ2 + B w2 (cf. (2.28)) and the foregoing expression, we obtain α(Ω) 2 First, we bound the last two terms in the same way as Lemma 2.3 (that is, using (2.26) and (2.30)): and with p, q ∈ [1, +∞) such that 1 p + 1 q = 1. We then take into consideration the further-regularity assumed in Section 2.4, and recall from the Sobolev embedding theorem that H ε (Ω) is continuously embedded into L 2p (Ω), with In this way, and (2.38) now yields (2.39) Therefore, putting (2.36), (2.37) and (2.39) together into (2.35), we obtain and since t 1 = S 1 (w 1 , φ 1 ) and u 1 = S 3 (w 1 , φ 1 ), the last inequality is exactly the estimate (2.34). We end the proof by noting that the foregoing expression defining C S confirms that this constant is independent of r.
We emphasize once again that, differently from [5], the present approach allows to handle both the 2D and 3D cases for the Boussinesq model with temperaturedependent viscosity. Indeed, notice how the foregoing Lemma shows the main difference with respect to [5]: the tensor-product term in (2.36) is no longer multiplied by an H 1 -term, thus avoiding the use of the injection H 1 (Ω) → L 8 (Ω) (not ensured for Ω ⊆ R 3 ) when splitting them as in [5, Eq. (3.54)], yielding a more robust formulation for the two and three-dimensional cases. On the other hand, the analogous result for S remains intact. Lemma 2.7. There exists a positive constant C S , independent of r, such that Proof. We refer to [20,Lemma 3.7] for full details. We just remark here that the resulting constant C S depends on the ellipticity constant of a in the kernel of b, and on the constant c(Ω) arising from the boundedness estimates of the functionals F w1,φ1−φ2 and F w1−w2,φ2 , namely (cf. [20, eq. (3.39)]): from which it is clear that C S does not depend on r.
As a consequence of the previous Lemmas, T is a Lipschitz continuous operator, as shown next.
Lemma 2.8. Let r ∈ (0, r 0 ), with r 0 given as in (2.29). Then, there exists a constant C T > 0 depending on r, such that Proof. The result comes from the definition of T (cf. (2.22)) and the estimates obtained in the previous two lemmas (cf. (2.34) and (2.40)) in an identical way to [20, Lemma 3.8] (see also [5,Lemma 3.10]). We omit further details.
In summary, from Lemmas 2.3 and 2.4, T : W → W is well-defined and does map the ball into itself (thanks to Lemma 2.5). Then, Lemma 2.8 shows that the operator T is Lipschitz continuous with a data-depending constant given by C T g ∞,Ω + u D 1/2+ε,Γ . When this Lipschitz constant is less than 1, the existence and uniqueness of a fixed point follows thanks to the Banach fixed-point theorem. In this regard, we stress that the dependence of C T on a given r does not affect the imposition of the aforementioned restriction since it is clearly achieved by sufficiently small data. By what has been explained in Section 2.2, this fact is equivalent to the well-posedness of the augmented mixed-primal formulation (2.11), and the result is stated as follows.
3. The Galerkin scheme. In this section we derive a Galerkin method for the augmented mixed-primal formulation (2.11). Let us consider T h , a regular triangulation ofΩ by triangles T (when n = 2) or tetrahedra T (when n = 3) of diameter h T and define the meshsize h := max{h T : T ∈ T h }. Then, consider arbitrary finite- Hence, according to (2.11), the Galerkin scheme reads: where the forms A ϕ h , B u h , a and b; and the functionals F ϕ h , F D , F u h ,ϕ h are defined by (2.12)- (2.19). Since the proof of well-posedness follows the steps of the previous section, and moreover, analogously to [5, Section 4], we only state the requirements to be imposed over the finite-dimensional subspaces and the main result, which is analogous to [5, Theorem 4.8].
3.1. Well-posedness of the Galerkin scheme. It can be seen that no restrictions have to be added to H t h , H σ h and H u h other than being finite-dimensional subspaces of the described spaces. However, for ellipticity purposes of a in the discrete kernel of the operator induced by b (according the the Babuška-Brezzi theory), the following two inf-sup conditions must be met: Then, denoting by W h the closed ball in H h := H u h × H ϕ h of radius r and center (0, 0), that is the main result of this section reads as follows.
Proof. We only mention that (3.1) is transformed into a fixed-point problem that is analyzed by means of the Brouwer fixed-point theorem in the convex and compact set W h (see [5,Theorem 4.8]).

3.2.
Specific choice of finite element subspaces. Given a set S ⊂ R := R n and an integer k ≥ 0, we define P k (S) as the space of polynomial functions on S of degree ≤ k, and for each T ∈ T h , we define the local Raviart-Thomas spaces of order k as RT k (T ) := P k (T ) + P k (T )x, where x is a generic vector in R. Hence, the strain rate, stress, velocity and temperature variables can be approximated using the following finite element subspaces: whereas for the normal component of the heat flux, we let { Γ 1 , Γ 2 , . . . , Γ m } be an independent triangulation of Γ (made of straight segments in R 2 or triangles in R 3 ), and define h := max j∈{1,...,m} | Γ j |. Then, with the same integer k ≥ 0 used in definitions (3.5)-(3.8), we approximate λ by piecewise polynomials of degree ≤ k over this new mesh, that is It can be seen that H ϕ h and H λ h satisfy the conditions (3.2) and (3.3) as long as h ≤ C 0 h for some C 0 > 0 (cf. [5,Section 4.3]). Other problems and corresponding references where some of or all the above finite element subspaces have been employed include, among others, coupled flow-transport [8,9], natural convection with phase-change [7], Navier-Stokes [15,16], and Stokes-Darcy (see, e.g. [29], where the above restriction between the mesh sizes h and h was discussed).

4.
A priori error analysis. Let ( t, (ϕ, λ)) ∈ H × H 1 (Ω) × H −1/2 (Γ) with (u, ϕ) ∈ W be the solution of (2.11), and ( t h , (ϕ h , λ h )) ∈ H h ×H ϕ h ×H λ h with (u h , ϕ h ) ∈ W h be a solution of the discrete problem (3.1), that is, and In what follows, we denote as usual First, the error estimate related to the variables of the momentum equation is obtained by means of the Strang Lemma, applied to the pair (4.1). We recall the Lemma, and its consequent result next.
is the ellipticity constant of A ϕ + B u (cf. (2.28)). Then, there holds Proof. From Lemma 2.3, we see that A ϕ +B u and A ϕ h +B u h are bilinear, bounded (both with constant A ϕ + B u , w.l.o.g. since it is independent of (u, ϕ)) and uniformly elliptic (both with constant α(Ω) 2 ). Also, F ϕ + F D and F ϕ h + F D are linear bounded functionals in H and H h , respectively. Hence, a straightforward application of Lemma 4.1 to the pair (4.1) yields Then, in order to estimate the last supremum in (4.4), we add and subtract suitable terms to write and so, using the boundedness of the bilinear forms A ϕ + B u and A ϕ h + B u h (cf. (2.27)), the estimate (2.8), the continuous embedding H 1 (Ω) → L n/ε (Ω) with constant C ε and the further-regularity assumption in a similar way to (2.39), we obtain This inequality, together with (4.5), back into (4.4), results in (4.3), therefore concluding the proof.
Hence, adding the estimates obtained in the previous two lemmas, we have a preliminary estimate for the total error: where Then, bounding the terms u 1,Ω , | ϕ | 1,Ω and u h 1,Ω using the continuous dependence results (2.41), (2.42) and (3.4), and the further-regularity assumption (2.33) to bound t ε,Ω , (4.6) becomes Therefore, using the continuous injection H 1/2+ε (Γ) → H 1/2 (Γ) with constant C i and denoting by C 5 := C 2 C S , C 6 := (C 1 + C 5 r + C 2 )C S , C 7 := C 6 r + C 3 C S (r)r + C 4 , we see from (4.7) that Then, there exists C > 0 depending only on parameters, data and other constants, all of them independent of h, such that Proof.
The assumption (4.9) allows us to subtract the total error term in the righthand side of (4.8), thus verifying the Céa's estimate with C = 2 max Finally, we state the rates of convergence of the Galerkin scheme (3.1) when the finite element subspaces (3.5)-(3.9) are used.

5.
Numerical results. We now present three numerical examples showing the performance of the augmented mixed-primal method (3.1) with the subspaces specified in Section 3.2. The computational implementation uses the finite element library FEniCS (cf. [6]), in particular the module multiphenics 1 that allows a straightforward description of the Lagrange multiplier and a block structure manipulation of the linearized systems arising at each Picard step. The solution of these linear systems is done employing the unsymmetric multi-frontal direct solver MUMPS (cf. [10]). The implemented iterative method mimics the fixed-point strategy from Section 2.2: it begins with an initial state at rest (u, ϕ) = (0, 0) and it stops whenever the relative error between two consecutive iterations of the solution vector measured in the 2 -norm (the sum of squared entries of a vector) is sufficiently small, i.e., where tol = 10 −6 is a specified tolerance. We also recall that the pressure is post-processed as which, as explained in [5,Section 5.2], converges to the exact pressure p at the same rate as the other variables (cf. Theorem 4.5). Similarly, as suggested by the second expression in (1.4), the vorticity can be postprocessed as which easily yields thus proving that it also converges at the same rate provided by Theorem 4.5. In this way, we define the error per variable Non-homogeneous terms will then appear in the right-hand sides of the momentum and energy equations. Nevertheless, well-posedness is still ensured, since the smoothness of the exact solution makes these terms immediately belong to L 2 (Ω), thus requiring only a minor modification in the functionals of the variational formulation. Concerning the stabilization parameters, these are taken as pointed out in Section 2.3, where the viscosity bounds are estimated as µ 1 = 0.5, µ 2 = 1.25, thus resulting in κ 1 = κ 2 = 0.32, κ 3 = 0.25 and κ 4 = 0.125. We also remark that the trace condition on the stress is enforced through penalization, here and also in the upcoming examples. In Figure 5.1 we show part of the obtained numerical solution with 214,788 DOF and a first order approximation, whereas in Table 5.1 we show the convergence history given the specified data and the finite element spaces from Section 3.2 with successive quasi-uniform mesh refinements. In both cases, it can be seen that the rates of convergence are the expected ones according to Theorem 4.5, that is, O(h) for the first case, and O(h 2 ) for the second one. We also observe that around eight iterations are needed to reach convergence of the Picard algorithm.   where r 1 and r 2 are positive parameters, and It is expected to find a counter-clockwise rotating vortex with center (x,ŷ), wherê x = 1 r 1 log e r1 + 1 2 ,ŷ = 1 r 2 log e r2 + 1 2 .
In particular, taking r 1 = r 2 = 4.5, the center of the vortex is expected to appear at the top-right corner of the cavity (that is, (x,ŷ) = (0.829, 0.829)). Then, considering  Table 5.2 we show the corresponding convergence history. As expected, when using the finite-element subspaces of Section 3.2 with k = 0 and k = 1, the rates of convergence are near the optimal linear and quadratic orders, respectively. Part of the solution is shown in Figure 5.2, where a second-order approximation has been used with 724,448 DOF. The second-order method is more convenient than the first-order scheme, in terms of CPU time used to reach the same precision. Notice that a relatively high refinement was required to capture the small features of the solution (e.g. pressure and pseudostress).     Table 5.3. Convergence history for Example 3, with a quasiuniform mesh refinement and approximations of first and second order.
Dirichlet datum on Γ. The exact temperature is uniformly bounded and it is also exploited as Dirichlet datum. In this configuration the viscosity bounds can be set as µ 1 = 0.5, µ 2 = 1 and the augmentation constants take the values κ 1 = κ 2 = 0.5, κ 3 = 0.25, and κ 4 = 0.125. The error history, associated with the schemes of order one and two, are performed using six steps of uniform mesh refinement applied to an initial structured tetrahedral mesh. On each level we compute approximate solutions, as well as errors and convergence rates defined as above. The boundary partition is considered conforming with the interior mesh, for sake of convenience and simplicity of the 3D mesh generation. Our findings are collected in Table 5.3, where errors and Picard iteration count are tabulated by number of degrees of freedom and meshsize. As in the 2D case, optimal error decay is observed for all individual errors (including the post-processed variables), and we also note that the errors e(σ), e(u) are dominant. One can also see that (perhaps assisted by the conformity between the interior and boundary meshes) the error associated with the boundary heat flux has a convergence systematically better than the optimal predicted by Theorem 4.5. Finally, we portray in Figure 5.3 a sample of the approximate solutions generated by the lowest-order mixed method on a relatively coarse mesh.

Example 4.
To close this section we conduct a benchmark test of a differentially heated, two-dimensional cavity of width 2 and height 1 (dimensionless units). The nonlinear viscosity used here is also explicitly dependent on depth, as in the models used in the context of mantle convection [12,40], µ(ϕ) = 2 exp(−9.7044ϕ + 4.1588(1 − y)), and the boundary treatment proceeds as follows. For the thermal energy conservation, the boundary is split into Γ D (top and bottom edges of the box) and Γ N (vertical walls) where temperature and heat flux are prescribed, respectively. The boundary temperature on Γ D is set to 0 on the top and 1 on the bottom edges, and it suffices with adequately defining ϕ D in the formulation. On Γ N we simply set zero thermal flux (the vertical walls are considered insulated). For the Navier-Stokes equations, and following [40], we impose free slip conditions defined as u · ν = σν · m = 0, where m denotes the tangential vector. These conditions imply that the Galerkin term for boundary velocity of Section 2 is now modified as κ 4 Γ (u · ν) (v · ν) = 0 ∀v ∈ H 1 (Ω), and we also require the additional condition κ 5 Γ (σν · m) (τ ν · m) = 0 ∀τ ∈ H(div; Ω), where κ 5 > 0 can be taken, e.g., equal to κ 4 . The coupled Boussinesq equations can be scaled so that the thermal conductivity is inverse to the adimensional Rayleigh number Ra=10 4 , and the buoyancy term is simply ϕg with g = (0, 1) t . In order to produce convection regimes we require extending the formulation to the transient case, adding to the thermal energy equation the term ∂ t ϕ and to the momentum equation the acceleration ∂ t u (which for mantle convection could eventually be dropped as it scales with the inverse of a very large Prandtl number). These time derivatives are approximated with the backward Euler scheme, using a constant timestep of ∆t = 0.1 and running the coupled Boussinesq equations until the final time t = 12. The initial velocity is zero and the initial temperature is ϕ(0) = 1 − y 2 + 0.01 cos(4πx) sin(πy). We employ a relatively coarse structured mesh with 32K triangles (but graded to be refined in the y−direction near the bottom and top boundaries), which results in a mixed-primal method (for the lowest-order case) having 208,986 DOF. For this test we observe that an average of five fixed-point iterations are required to reach the fixed tolerance of 10 −7 . The Rayleigh number used here was sufficiently large to generate thermal blob-shaped forms (see the plumes in the temperature plot of Figure 5.4). The obtained convection modes are consistent with the flow patterns observed in [40] (see also [31,33]). They indicate that the velocity is driven by the temperature difference.