Divergence-conforming methods for transient doubly-diffusive flows: A priori and a posteriori error analysis

The analysis of the double-diffusion model and $\mathbf{H}(\mathrm{div})$-conforming method introduced in [B\"urger, M\'endez, Ruiz-Baier, SINUM (2019), 57:1318--1343] is extended to the time-dependent case. In addition, the efficiency and reliability analysis of residual-based {\it a posteriori} error estimators for the steady, semi-discrete, and fully discrete problems is established. The resulting methods are applied to simulate the sedimentation of small particles in salinity-driven flows. The method consists of Brezzi-Douglas-Marini approximations for velocity and compatible piecewise discontinuous pressures, whereas Lagrangian elements are used for concentration and salinity distribution. Numerical tests confirm the properties of the proposed family of schemes and of the adaptive strategy guided by the {\it a posteriori} error indicators.


Introduction and problem formulation
1.1. Scope. A number of physical problems of relevance in industrial applications involve coupled incompressible flow and double-diffusion transport. We are interested in numerical schemes for the approximation of a class of coupled equations that arise as models of sedimentation of small particles under the effect of salinity of the ambient fluid. The governing model can be stated as follows (cf., e.g., [19,33]): posed on a spatial domain Ω ⊂ R d , d = 2 or d = 3, where t ∈ (0, t end ] is time, u is the fluid velocity, ν is the concentration-dependent viscosity, ρ m is the mean density of the fluid, p is the fluid pressure, ρ is density, g is the gravity acceleration, s is the salinity concentration, c is the concentration of solid particles, Sc = ν ref /κ s is the Schmidt number, where κ s is the diffusivity of salinity, ν ref is a reference viscosity in the absence of solid particles, and τ = κ s /κ c is the inverse of the diffusivity ratio, where κ c is the diffusivity of solid particles, and e z is the upward-pointing unit vector. We relate the densities through a linearised equation of state ρ = ρ m (αs + βc).
Again as in [19,33], the solid particles are assumed mono-sized with radius r, and settle at dimensionless veloc- is the Stokes velocity (settling velocity of a single particle in an unbounded fluid). The coupling mechanisms between flow and transport are only due to advection for concentration and salinity (where the advecting velocity for concentration, u − v p e z , is also divergence-free), and through the concentration-dependent viscosity. Further details are provided in later parts of the paper.
To put the paper further into the proper perspective, we mention that there exists an abundant body of literature devoted to constructing accurate finite element and related schemes for doubly-diffusive flows. Some recent contributions include variational multiscale stabilised schemes, least-squares methods, divergence-conforming mixed methods, volume-averaging discretisations, spectral elements, vorticity-based finite element formulations, and similar methods applied to, e.g., flows with heat and mass transport [1], reactive Boussinesq flows [2], nonlinear advection-reaction-diffusion in the context of bioconvective flows [10,31], cross-diffusion and boundary layer effects in doubly-diffusive Navier-Stokes-Brinkman equations [18,20] and in Darcy-Brinkman equations [39], or phase change models [21,38,40,41]; where the list is far from exhaustive.
The solvability analysis for the continuous and discrete problems usually follows energy and fixed-point schemes as done for classical Boussinesq equations, and this is also the approach we follow here. The discretisation in space uses an interior penalty divergence-conforming method for the flow equations (in this case, Brezzi-Douglas-Marini (BDM) elements of degree k ≥ 1 for the velocity and discontinuous elements of degree k − 1 for the pressure following [11,29]), combined with Lagrangian elements for the diffusive quantities, and the development stands as a natural extension of the formulation in [18] to the transient case. As such, it also features exactly divergencefree velocity approximations ensuring local conservativity and energy stability, and the error estimates of velocity are pressure-robust. The chosen time discretisation is the backward differentiation formula of degree 2 (BDF2), which for k = 2 gives a method of order 2 in space and time. Existence of discrete solutions is established by the Brouwer fixed-point theory similarly as in [18], and the error analysis in the semi-discrete and fully-discrete settings is adapted from the theory of [5] for the Boussinesq equations.
In many applications where double-diffusion effects occur, complicated flow patterns exist in zones far from boundary layers and sufficiently refined meshes are needed essentially in the whole spatial domain. However, for salinity-driven settling of solid particles that result in mathematical models such as (1.1), many of the flow features are clustered near zones of high-gradients of concentration, which is where the typical plumes are observed [19,31]. This motivates the use of adaptive mesh refinement guided by a posteriori error indicators. For instance, in the context of phase change models there are some results based on error-related metric change [21,34] and on goal-oriented adaptivity [41]. Regarding the design and rigorous analysis of residual-based a posteriori error estimators for flow-transport couplings, the literature is predominantly focused on the stationary case (see, e.g., [3,6,7,9,13,22,37,42] and the references therein). Only a few results are available for the time-dependent regime, from which we mention the adaptive mixed method for Richards equation in porous media [15], the remeshing scheme based on goal-oriented adaptivity for solidification problems advanced in [14], the collection of adaptive schemes for reactive flow discussed in [16] and for heat transfer in [30]. However, none of these theoretical frameworks is directly applicable to (1.1) using divergence-conforming approximations.
The a posteriori error analysis we advance here is of residual type, and its analysis uses ideas from the abstract results related to spatial estimators for discontinuous Galerkin schemes applied to parabolic problems in [24]. The approach hinges on a decomposition of the discrete solution into a conforming and a non-conforming contribution, along with a reconstruction technique (see also [32]). This has also been exploited for the construction of a posteriori estimators of time-dependent Stokes and Navier-Stokes equations [12,43]. Our a posteriori error analysis is divided into three parts. In first part, we present the error estimator for the steady coupled problem. In second part, we extend the a posteriori error estimation to the semi-discrete method, and finally we present the a posteriori error estimator for the unsteady coupled problem. We restrict that part of the analysis to the simpler backward Euler method. To the best of our knowledge, the a posteriori error estimation advanced in this paper is the first comprehensive study targeted for transient doubly-diffusive flows.
The remainder of this paper is organised as follows. In what is left of this section we outline the weak formulation of (1.1) and state the stability the continuous problem. In Section 2 we introduce the divergenceconforming method in fully discrete form, show existence of discrete solutions using fixed-point arguments, and rigorously establish a priori error estimates. Section 3 is devoted to the construction and analysis of efficiency and reliability for a residual-based a posteriori error estimator tailored for the stationary problem. In turn, these upper and lower bounds are used to establish properties of a second family of estimators for the transient case, and addressed in Sections 4 and 5. In Section 6 we collect numerical tests that verify the theoretical convergence rates predicted by the a priori error analysis, confirm the robustness of the proposed a posteriori error estimators, and illustrate the advantages of adaptive methods in the simulation of doubly-diffusive flows.
A problem similar to (1.2) has been studied in [4]. Assuming that F belongs to L 2 (0, t end ; H −1 (Ω)), that the initial velocity u 0 belongs to L 2 (Ω) and the initial data for the coupled species (s, c in the context of our problem) belongs to L 2 (Ω), the authors showed existence of a solution by using the Galerkin method and applying the Cauchy-Lipschitz theorem, and proved its uniqueness in two dimensions. Such analysis can be applied to (1.2) by noting that F is a Lipschitz-continuous function, and assuming the initial data belongs to appropriate spaces. This is, however, not the focus of the paper.

Finite element discretisation and a priori error bounds
We discretise space by a family of regular partitions, denoted T h , of Ω ⊂ R d into simplices K (triangles in 2D or tetrahedra in 3D) of diameter h K . We label by K − and K + the elements adjacent to an edge, while h e stands for the maximum diameter of the edge. If v and w are smooth vector and scalar fields defined on T h , then (v ± , w ± ) denote the traces of (v, w) on e that are the extensions from the interior of K + and K − , respectively. Let n e denote the outward unit normal vector to e on K. The average {{·}} and jump · operators are defined as whereas for boundary jumps and averages we adopt the conventions {{v}} = v = v, and {{w}} = w = w. In addition, ∇ h will denote the broken gradient operator.
For k ≥ 1 and a mesh T h on Ω, let us consider the discrete spaces (see e.g. [17]) [29]). Here P k (K) denotes the local space spanned by polynomials of degree up to k and V h is the space of divergence-conforming BDM elements. We then state the following semi-discrete Galerkin formulation for problem (1.2) The discrete versions of the trilinear forms a h 1 (·; ·, ·) and c h 1 (·; ·, ·) are defined using a symmetric interior penalty and an upwind approach, respectively (see, e.g., [11,29]): and u e h is the trace of u h taken from within the exterior of K.
We partition the interval [0, t end ] into N subintervals [t n−1 , t n ] of length ∆t. We use the implicit BDF2 scheme where all first-order time derivatives are approximated using the centred operator (similarly for ∂ t c) and for the first time step a first-order backward Euler method is used from t 0 to t 1 , starting from the interpolates u 0 h , s 0 h of the initial data. In what follows, we define the difference operator Dy n+1 := 3y n+1 − 4y n + y n−1 for any quantity indexed by the time step n. For instance, (2.3) can be written as The resulting set of nonlinear equations is solved by an iterative Newton-Raphson method with exact Jacobian. Hence the complete discrete system is given by For the subsequent analysis, we introduce for r ≥ 0 the broken H r (T h ) space as well as the following mesh-dependent broken norms where the stronger norm · 2,T h is used to show continuity. This norm is equivalent to · 1,T h on H 1 (T h ) (see [11]). Finally, adapting the argument used in [28,Proposition 4.5] we have the discrete Sobolev embedding: for r = 2, 4 there exists a constant C emb > 0 such that Using these norms, we can establish continuity of the trilinear and bilinear forms involved. The proof follows from [11,Section 4].
Lemma 2.1. The following properties hold: and for all u, v, w ∈ H 1 (T h ) and ϕ, ψ ∈ H 1 (Ω), there holds where the constantC Lip > 0 is independent of h (cf. [18]). Note that while the coercivity of the form a 2 (·, ·) in the discrete setting is readily implied by (1.5b), there also holds (cf. [29,Lemma 3 provided that the stabilisation parameter a 0 > 0 in (2.2) is sufficiently large and independent of the mesh size.
Let w ∈ H 0 (div 0 ; Ω), then due to the skew-symmetric form of the operators c h 1 and c 2 , and the positivity of the non-linear upwind term of c h 1 , we can write as well as the following relation (which is based on (2.5) and follows by the same method as in [28]): For any We also have Finally, we recall from [29] the following discrete inf-sup condition for b(·, ·), whereβ is independent of h: We will use the following algebraic relation: for any real numbers a n+1 , a n , a n−1 and defining Λa n := a n+1 − 2a n + a n−1 , we have 2(3a n+1 − 4a n + a n−1 , a n ) = |a n+1 | 2 + |2a n+1 − a n | 2 + |Λa n | 2 − |a n | 2 − |2a n − a n−1 | 2 .
(2.12)  We proceed in the same way taking ψ h = 4c n+1 h in the fourth equation of (2.4), to get the third result. We get the first result by substituting bounds for c h and s h into (2.14).
Theorem 2.2 (Existence of discrete solutions). The problem (2.4) admits at least one solution Proof of Theorem 2.2. To simplify the proof we introduce the following constants: ,Ω . We proceed by induction on n ≥ 2. We define the mapping

Note this map is well-defined and continuous on
and employ (2.9a), (2.9b) and (2.8), we obtain Next, using (2.13) we deduce that Then, setting Hence, the right-hand side is nonnegative on a sphere of radius r := C r /C R . Consequently, by Theorem 2. Let us denote by I h : C(Ω) → M h the nodal interpolator with respect to a unisolvent set of Lagrangian interpolation nodes associated with M h . Π h u will denote the BDM projection of u, and L h p is the L 2 -projection of p onto Q h . Under usual assumptions, the following approximation properties hold (see [29]): The following development follows the structure adopted in [5].
Lemma 2.2. Assume that u ∈ H 2 (Ω), p ∈ L 2 (Ω) and s, c ∈ H 1 (Ω). Then we have Proof. Since we assume u ∈ H 2 (Ω), integration by parts yields the required result. See also [11]. The third and fourth equations are a straightforward result from the continuous form.
Now we decompose the errors as follows: Assuming that u 0 h = Π h u(0), s 0 h = I h s(0) and c 0 h = I h c(0), we will use also the notation E n u = (u(t n )−Π h u(t n )) and ξ n u = (Π h u(t n ) − u n h ), and similar notation for other variables. Since for the first time iteration of system (2.1) we adopt a backward Euler scheme, an error estimate is required for this step.
and also that u L ∞ (0,t end ;W 1,∞ (Ω)) < M for a sufficiently small M > 0 (a precise condition can be found in Theorem 2.5). Then there exist positive constants C 1 u , C 1 s , C 1 c , independently of h and ∆t, such that Proof. As in the continuous case, we define the discrete kernel of the bilinear form b(·, ·) as and relying on the inf-sup condition (2.11), we can continue with an equivalent discrete problem without pressure.
Taking into account the assumed regularity for u, we have for all x, a γ(x) ∈ (0, 1) such that then u satisfies the following error equation , which results after choosing ξ 1 u as test function in the first equation of the reduced form of Lemma 2.2 and system (2.1), performing an Euler scheme step, subtracting both equations, and adding ±a h 1 (c 1 h ; Π h u(∆t), ξ 1 u ). Now, invoking the approximation estimates (2.15), Young's inequality, and the stability properties, we get (2.17) Next, we choose ξ 1 s as test function in (2.16c) and system (2.1); we follow the same steps as before, adding to the sum of both equations the term ±a 2 (I h s 1 , ξ 1 c ) to obtain 1 4 ξ 1 In the same way, choosing ξ 1 c as test function in (2.16d) and in (2.1) we obtain Now, from (2.17) we deduce that and u L ∞ (0,t end ;H 1 (Ω)) < M for a sufficiently small constant M > 0. Then there exist positive constants C, γ 1 ≥ 0 independent of h and ∆t such that for all m + 1 ≤ N , Proof. We appeal to the reduced form of the problem again, taking u h ∈ X h and u ∈ X, then we choose as test function v h = ξ n+1 u in the first equation of (2.4) and insert the terms We consider (2.16a) (see Lemma 2.2) at t = t n+1 and v = ξ n+1 u . Inserting the term ±(Du(t n+1 ), ξ n+1 u ) Ω /(2∆t), we readily deduce the identity We can then subtract (2.23) from (2.22) and multiply both sides by 4∆t, yielding I 1 + I 2 + · · · + I 7 = 0, with Du(t n+1 ), ξ n+1 u Ω , .
For the first term, using (2.12) we can assert that ,Ω . Using the ellipticity stated in (2.8), we readily get By using Taylor's formula with integral remainder we have then by combining Cauchy-Schwarz and Young's inequality, we obtain the bound Now we insert ±4∆tE u (t n+1 ) into the fourth term, which leads to Proceeding as before and using (2.15) on the first term of I 4 , we get . Then by (2.15), (2.7), Lemma 2.1 and Young's inequality we immediately have Adding and subtracting suitable terms within I 6 yields where we defineĨ The bound (2.10) and (2.15) imply that where C * is a positive constant coming from (2.15). We also have Hence, by choosing ε i = 3α a /5 for i = 1, . . . , 10, collecting the above estimates, and summing over 1 ≤ n ≤ m for all m + 1 ≤ N , we get where Finally, Theorem 2.4 yields the desired result.
Theorem 2.6. Let (u, p, s, c) be the solution of (1.2) under the assumptions of Section 1.4, and (u h , p h , s h , c h ) be the solution of (2.4). If Proof. Proceeding similarly as for Theorem 2.5, we choose as test function ϕ h = ξ n+1 s in the second equation of (2.4) and insert suitable additional terms to obtain the following identity (analogous to (2.22 From (2.16b), focusing on t = t n+1 , using ϕ = ξ n+1 s and proceeding as in the derivation of (2.23), we obtain Next we subtract (2.25) from (2.24) and multiply both sides by 4∆t. This yieldsÎ 1 + · · · +Î 6 = 0, wherê Inserting ±4∆tE s (t n+1 ) intoÎ 4 and using (2.15) leads to the bound Employing again (2.15) in combination with (1.3a) we have In order to derive a bound forÎ 6 we proceed as for the bound on I 7 in the proof of Theorem 2.5; namely adding and subtracting suitable terms in the definition ofÎ 6 , definingĨ 6 in this case bŷ and applying (2.9b), (2.6a), (2.15) and Young's Inequality to the result, we get In this manner, and after choosing ε i = 6α a /7 for i = 1, . . . , 7, we can collect the above estimates and sum over 1 ≤ n ≤ m, for all m + 1 ≤ N , to get This concludes the proof.  Proof. It follows along the same lines of the proof of Theorem 2.6, with constant γ 3 given by Proof. Owing to the inf-sup condition (2.11), there exists From (2.4) and Lemma 2.2, proceeding as in the proof of Theorem 2.5, we obtain Summing over 1 ≤ n ≤ m for all m + 1 ≤ N and substituting back into (2.27), we obtain , and the desired result readily follows from Theorem 2.8.
We next proceed to derive and analyse a posteriori error estimators. We split the presentation into three cases of increasing complexity, starting with an estimator focusing on the steady coupled problem.
Next, we decompose the H(div)-conforming velocity approximation uniquely into Proof. It follows straightforwardly from the decomposition u h = u c h + u r h and from the edge residual. Lemma 3.2. If u 1,∞ < M , s ∞ < M and c ∞ < M , then the following estimate holds: Proof. Using u h = u c h + u r h and the triangle inequality imply

Owing to the relation
, while using the properties and we readily see that after stating the discrete problem as and employing (3.4), the sought results follow. (3.5)

Proof. Using integration by parts
gives where we define the terms Applying the Cauchy-Schwarz inequality to T 1 implies Next, we rewrite T 2 in terms of a sum over interior edges and apply again the Cauchy-Schwarz inequality. Then Then, owing to the Cauchy-Schwarz inequality, it follows that Proceeding similarly, we may establish the following bounds for T 4 and T 5 : C ∇φ 0,Ω , where C > 0 is a constant independent of h.
Proof. It suffices to apply Lemmas 3.2 and 3.3.

3.2.
Efficiency. For each K ∈ T h , we can define the standard polynomial bubble function b K . Then, for any polynomial function v on K, the following results hold: where C is a positive constant independent of K and v.
Lemma 3.4. The following estimates hold, where C is a positive constant: Moreover, it also follows that Then, using (3.10), we have Using the Cauchy-Schwarz inequality and (3.10) we obtain , and combining these bounds leads to the first stated result. The other two bounds follow similarly.
Let e denote an interior edge that is shared by two elements K and K . Let ω e be the patch which is the union of K and K . Next, we define the edge bubble function ζ e on e with the property that it is positive in the interior of the patch ω e and zero on the boundary of the patch. From [36], the following results hold: Moreover, we also have Proof. Let e be an interior edge and let us define a rescaling of the edge bubble function in the form Using (3.11) gives (3.12) Using integration by parts on each element of patch ω e implies Note that (u, p, s, c) solves the underlying problem, so we then have Next, applying the Cauchy-Schwarz inequality together with Lemma 3.4 and (3.11) gives Combining the bounds of T 1 , T 2 and T 3 with (3.12) and (3.13) implies the first stated result. Similarly, we can prove the other two bounds. Let Ψ be defined as in (3.3). Then there exists a constant C > 0 that is independent of h such that Proof. Combining Lemmas 3.4 and 3.5 implies the stated result.

A posteriori error bound for the semidiscrete method
For each t ∈ (0, T ], let us consider the problem: find (ũ,p,c,s) Also, for each t ∈ (0, T ], we write the discrete weak formulation: where (4.1) remains in effect.
Next we introduce the error indicator Θ as where whereas Ψ is the global a posteriori error estimator for the steady problem with element and edge residual contributions defined in (3). In this case we now replace f and f 1 , f 2 by (4.1).  On the other hand, integrating by parts in Lemma 4.1 yields − a 1 (c,ũ, v) − c 1 (e u ;ũ, v), ∀v ∈ H 1 0 (Ω) (∂ t e s , φ) = −a 2 (ρ s , φ) − c 2 (u; ρ s , φ) − c 2 (e u ;s, φ), ∀φ ∈ H 1 0 (Ω), We apply Young's inequality and the definition of the dual norm. Then, we integrate in time the resulting expression. Finally, the second result is a consequence of Theorem 3.2 and (4.4).

A posteriori error analysis for the fully discrete method
In this section, we develop an a posteriori error estimator for the fully discrete problem and focus the presentation on the simpler case of a time discretisation by the backward Euler method. For each time step k (1 ≤ k ≤ N ), we define the (global in space) time indicator Ξ k as Next we define the accumulated time and spatial error indicators as where the terms Υ 2 k are constructed with the a posteriori error estimator contributions defined as in the steady case (3), but at a given time step k. That is, For each time step k, we can split again the H(div)-conforming discrete solution u k h into a conforming part u k hc and a non-conforming part u k hr such that u k h = u k h,c + u k h,r . For each t ∈ (t k−1 , t k ], we introduce a linear interpolant u h (t) in terms of t as where {l k , l k+1 } is the standard linear interpolation basis defined on [t k , t k+1 ]. Similarly, we may introduce u h,c (t) and u h,r (t). Then, setting e uc = u − u h,c , we have e u = u − u h = e uc − u h,r . For t ∈ (t k−1 , t k ), we define and for all t ∈ (t k−1 , t k ), we consider the problem of finding (ũ k ,p k ,s k , Lemma 5.1. The following estimates hold Proof. Combining (1.2) and (5.2) implies Choosing v = e uc τ , q = p − p k , φ = e s τ , ψ = e c τ and then combining the first two equations, we have (∂ t e uc τ , e uc τ ) + a 1 (c, e uc τ , e uc τ ) + c 1 (u, e uc τ , e uc τ ) − (∂ tũh,r , e uc τ ) + a 1 (c − c h ,ũ, e uc τ ) + c 1 (u − u h , u, e uc τ ) + a 1 (c,ũ h,c −ũ k , e uc τ ) + c 1 (u,ũ h,c −ũ k , e uc τ ) = 0, (∂ t e s τ , e s τ ) + a 2 (e s τ , e s τ ) + c 2 (u; e s τ , e s τ ) − c 2 (u − u h ;s k , e s τ ) + a 2 (s h −s k , e s τ ) − c 2 (u; s h −s k , e s τ ) = 0, where In light of the definition of uτ , sτ and cτ , we get Then we can apply triangle inequality, which gives Finally, applying integration by parts in (5.3a) yields Next we apply Young's inequality and the definition of the dual norm. Then, we integrate the whole expression in time between t k−1 and t k for each k = 1, 2, . . . , n and sum the expression for each k. Finally, we use (5.4), (5.5), (5.6) and (5.7) to establish the second stated result.
Proof. Using u k h = u k h,c + u k h,r together with the identity in [24, (5.59)-(5.60)] that in our context reads , we can invoke Lemma 5.1 and reuse the strategy applied in Theorem 4.1 to complete the proof.
As u is prescribed everywhere on ∂Ω, for sake of uniqueness we impose p ∈ L 2 0 (Ω) through a real Lagrange multiplier approach. To verify the a priori error estimates, we introduce the discrete norms where e,ẽ denote errors generated on two consecutive pairs of mesh size and time step (h, ∆t), and (h,∆t), respectively. Choosing ∆t = √ 2h and using scheme (2.4), the results in Table 6.1 confirm that the rates of convergence are optimal, coinciding with the theoretical bounds anticipated in Theorem 2.8. Approximate velocity magnitude (after 3 refinement steps), pressure (after 4 refinement steps), concentration s (after 5 refinement steps), and distribution of c after 6 steps of adaptive refinement.
6.2. Example 2: adaptive mesh refinement for the stationary problem. The classical Dörfler strategy [23] is employed for the adaptive algorithm based on the steps of solving, estimating, marking, and refining. Estimation is performed by computing the error indicators and using them to select/mark elements that contribute the most to the error [30]. The marking is done following the bulk criterion of selecting sufficiently many elements so that they represent a given fraction of the total estimated error. That is, one refines all elements K ∈ T h for which where 0 < γ ratio < 1 is a user-defined constant (that we tune in order to generate a similar number of degrees of freedom, or comparable errors, as those obtained under uniform refinement). And then the algorithm aims for equidistribution of the local error indicator in the updated mesh.
In the adaptive case, instead of (6.1) the convergence rates (for the spatial errors) are computed as where DoF and DoF are the number of degrees of freedom associated with each refinement level. The robustness of the global estimators is measured using the effectivity index (ratio between the total error and the indicator) .
These solutions exhibit a generic singularity towards the reentrant corner, and therefore one expects that the error decay is suboptimal when applying uniform mesh refinement. After solving the coupled stationary problem on sequences of uniformly and adaptively refined meshes, the aforementioned behaviour is indeed observed in Table 6.2, where the first part of the table shows deterioration of the convergence due to the high gradients of the exact solutions on the non-convex domain. The results shown in the bottom block of the table confirm that as more degrees of freedom are added, a restored error reduction rate is observed due to adaptive mesh refinement guided by the a posteriori error estimator Ψ. The second-last column of the table also indicates that the effectivity index oscillates under uniform refinement, while it is much more steady in the adaptive case. We tabulate as well the Newton-Raphson iteration count (needed to reach the relative residual tolerance of 1e-6), and this number is also systematically smaller for the adaptive case (about four steps in all instances) than for the uniform refinement case (up to seven nonlinear steps for certain refinement levels). As an example we plot  c(x, y, t) = 1 2 (1 + cos(π/4(xy))) exp(−t), s(x, y, t) = 1 2 (1 + sin(π/2(xy))) exp(−t).
Cumulative errors up to t final are computed as , and the resulting error history, after six steps of uniform mesh refinement, is collected in Table 6.3. To be consistent with the development in Section 5, the numerical verification in this set of tests was carried out using a backward Euler time discretisation. The a posteriori error estimator (5.1) is computed and the effectivity index is also tabulated, showing that the estimator is robust (and confirming the theoretical reliability bound as well as giving an heuristic indication of its efficiency). Note that in this case, since the mesh refinement is uniform, the auxiliary interpolation of the solutions at the last time step on the current mesh is not necessary. The average number of Newton-Raphson iterations required for convergence was 3.2.
6.4. Example 4: simulation of salinity-driven flow instabilities. To illustrate the behaviour of the model and the proposed method, we simulate a salinity-driven flow problem. Similar examples are found in [19] where direct numerical simulations (DNS) are applied to the version of (1.1) that has constant viscosity.
The configuration of layering in sedimentation is taken from [33]. We consider a rectangular domain with L x = 40, L y = 300 and an initial solid-particle concentration profile that is periodic in the horizontal direction, and periodic with Gaussian noise in the vertical direction: s(x, y, z, 0) = A 0 exp(−z 2 /σ 2 ) + A 1 sin(x),  with initial amplitudes A 0 , A 1 and width σ (see Figure 6.2). For the velocity field, we use non-slip boundary condition in all four walls and we choose ∆t = 0.1. As discussed in [33], simulations at low density ratios are extremely costly because of the large Reynolds numbers of fingering convection. In consequence we choose an initial density ratio R 0 = αs 0z /βc 0,z ≈ 4, and we carry out the simulations on tall, thin domains. Apart from the specifications above, the remaining constant parameters needed in the model take the following values A 0 = 2.86, A 1 = 0.5, σ = 0.35, ν = 1 × 10 −3 kg/m 3 , g = 9.8m/s 2 , Sc = 7.0, τ = 25, v p = 0.04m/s, α = −2.0, β = 0.5.
According to [33] a linear fingering instability occurs provided 1 < R 0 < τ , hence the instability shown in Figure 6.2 is expected. 6.5. Example 5: adaptive simulation of exothermic flows. To conclude this section, and to include an illustrative simulation exemplifying that the H(div)-conforming scheme along with the a posteriori error estimator perform well for an applicative problem, we address the computation of exothermic flows that develop fingering instabilities. The problem configuration is adapted as a simplification of the problem solved in [31], where the fields c, s represent solutal concentration and temperature, respectively. The model assumes an additional drag term due to porosity so that the momentum equation is of Navier-Stokes-Brinkman type. The domain is the rectangular region Ω = (0, L) × (0, H), and the initial solutal and temperature profiles are imposed as Boundary conditions are of mixed type for solutal and temperature. Both fields are prescribed to 0 and 1 on the bottom and top of the domain, respectively; while on the vertical walls we impose zero-flux boundary conditions. The velocity is of slip type on the whole boundary, and therefore a zero-mean condition for the pressure is considered using a real Lagrange multiplier. The solution algorithm, differently from the previous tests, is based on an inner fixed-point iteration between an Oseen and a transport system, rather than an exact Newton-Raphson method. An initial coarse mesh of 5300 elements is constructed, and an adaptive mesh refinement (only one iteration) guided by the estimator (5.1) is applied at the end of each time step. Figure 6.3 shows snapshots of adapted meshes at different times, and also samples of solute concentration, temperature distribution, velocity, and pressure at the final time.