COUPLED MIXED FINITE ELEMENT AND FINITE VOLUME METHODS FOR A SOLID VELOCITY-BASED MODEL OF MULTIDIMENSIONAL SEDIMENTATION ⋆

. In this paper we introduce and analyze a model of sedimentation based on a solid velocity formulation. A particular feature of the governing equations is given by the fact that the velocity field is non-divergence free. We introduce extra variables such as the pseudostress tensor relating the velocity gradient with the pressure, thus leading to a mixed variational formulation consisting of two systems of equations coupled through their source terms. A result of existence and uniqueness of solutions is shown by means of a fixed-point strategy and the help of the Babuˇska–Brezzi theory and Banach theorem. Additionally, we employ suitable finite dimensional subspaces to approximate both systems of equations via associated mixed finite element methods. The well-posedness of the resulting coupled scheme is also treated via a fixed-point approach, and hence the discrete version of the existence and uniqueness result is derived analogously to the continuous case. The above is then combined with a finite volume method for the transport equation. Finally, several numerical results illustrating the performance of the proposed model and the full numerical scheme, and confirming the theoretical rates of convergence, are presented.

With the help of dimensional analysis many models are reduced to a scalar transport-flow equation for the local concentration of solid particles coupled with a Stokes-type system for the velocity-pressure pair.
Inspired by the work done by Gustavsson and Oppelstrup [46], we study a model of multidimensional sedimentation based on the solid phase velocity field , pressure and volume fraction of the solid particles . The model equations, obtained by combining the ideas of [16,46], are essentially a variant of the model presented in [46], where the velocity field has the peculiarity of not being divergence-free. The system is strongly coupled containing nonlinear terms, and therefore the construction of numerical schemes for its approximation are involved and need a special treatment. For an introduction to the theory of sedimentation based on the postulates of continuum mechanics for mixtures of two superimposed continuous media we refer to [24,35,63]. One-dimensional PDE models of sedimentation have been studied widely in the literature, from the development of numerical approximations for batch and continuous settling to the study of steady state solutions, see e.g., [15, 18, 19, 31-33, 37, 41]. The special case of quasi-one-dimensional models where the influence of the domain geometry is reflected in coefficient terms of a one-dimensional conservation law have been studied in [8,17,22,23,25].
Several mathematical models and numerical approximations of multidimensional sedimentation have been reported in the literature. Bürger et al. [16] presented a complete derivation of a multidimensional model of sedimentation starting from the balance equations (mass and momentum) for the solid and fluid phases. After a dimensional analysis the model is reduced to a Stokes-type system coupled to a strongly degenerate convection-diffusion equation, consisting of three unknowns and whose velocity field is the so called "volume average velocity of the mixture" denoted by . Borhan and Acrivos [14] introduce a two-dimensional model for the batch sedimentation in inclined channels, with a single momentum equation for the mixture which varies in its definition depending on the fluid-suspension interface. In [50], a detailed discussion about the viscous stresstensors of the solid phase and mixture based on the rheological complexity of colloidal suspensions is provided. They also show numerical simulations of their approach in cylindrical coordinates for the case of a "deep-cone" thickener. Rao et al. [55] presented a two-phase model based on the so called "mass-average velocity" (denoted by ) for the simulation of batch sedimentation and shear between concentric rotating cylinders. The velocity field does not satisfy the condition of divergence free, and the continuity equation is no longer solenoidal since the density varies with the local volume fraction. The batch model derived by Gustavsson and Oppelstrup [46] is based on the solid phase velocity field and contemplates the possibility of a horizontal movement of the bottom of the vessel by having Dirichlet boundary conditions on . One of the main features of this model that differs from typical -based ones (see [16,21]) is that like the model in [55], the divergence of the velocity field is non-zero but varies with the pressure and local volume fraction. A finite differences scheme is introduced in [46] to simulate the model equations. Up to the authors knowledge, the well-posedness of the decoupled system for the velocity-pressure pair has not been studied yet. Regarding well-posedness of the fully coupled system, we refer to the analysis made for a linearized version of the model in [45].
A variety of advanced numerical schemes have been used to compute numerical approximations for multidimensional models of sedimentation [20,21,47,53,55,62]. In [47], the authors develop a finite element numerical scheme for a -based model coupled with the " -equations" for turbulent flows. They also treat the time approximations with a fractional-step scheme in order to have second order accuracy and preserve physical oscillations. Bürger et al. [20], introduced a multiresolution finite volume scheme for the two-dimensional batchsedimentation model given in [16], and presented numerical examples of the process in channels with inclined walls. In [21], the authors worked on an advanced numerical method combining a stabilized finite element scheme with finite volumes in staggered grids for the continuous sedimentation of the model [16]. Numerical simulations of the continuous sedimentation applied to wastewater treatment plants are also included. Another paper dedicated to the approximation of such models using a combined finite element-finite volume schemes can be found in [53]. In [62], a CFD approach for a -based model coupled to -equations is considered. On the other hand, and regarding the coupled flow and transport problem driven by a scalar non-linear reaction-diffusion equation interacting with the Stokes or Brinkman equations, which models the transport of species within a viscous flow, we stress that diverse combinations of primal and mixed finite element methods have been developed in recent years for its numerical solution (see, e.g., [3][4][5][6]11,12,27]). Mathematical modeling based on CFD simulations of two-dimensional sedimentation can be found in [26,29,36,48,49]. Experiments of sedimentation in vessels with inclined walls can be found in [57,60], and studies on the Boycott effect produced in these vessels were reported in [58,59]. Sedimentation models coupled to an energy balance and a temperature equation are studiend in [54,56]. Applications of sedimentation models can be found in mineral processing [30,34,38,39], wastewater treatment [1,40,62], petroleum industry [51,52], magma chambers [13,61] and oil-in-water emulsions [28].
According to the previous discussion, the first purpose of this paper is to analyze the existence and uniqueness of weak solutions for the decoupled velocity-pressure system. Secondly, we seek to propose reliable finite element methods for its numerical solution, and finally we aim to approximate the transport equation by a conservative finite volume scheme. The paper is organized as follows. The rest of this section collects some useful notation to be employed along the paper. Then, in Section 2 we derive the model equations. Next, in Section 3 we introduce some auxiliary unknowns to derive a mixed variational formulation for the aforementioned system, and employ a fixed-point strategy along with the Babuška-Brezzi theory and Banach theorem to address its solvability. In Section 4, we consider suitable finite dimensional subspaces to define the associated mixed finite element scheme, and adopt the discrete analogue of the analysis from Section 3 to prove its well-posedness. A priori error estimates and corresponding rates of convergence are also provided here. Then, in Section 5 we propose a finite volume discretization for the transport equation and couple the resulting procedure with the velocity-pressure mixed finite element formulation from Section 4. Finally, illustrative numerical results are shown in Section 6 and concluding remarks are presented in Section 7.

Further notation
In this section we recall some standard notation used in the rest of the paper. For any vector field = ( ) =1, in R we let ∇ and div( ) be the gradient and divergence of , respectively. In addition, for any tensor = ( ) , =1, we let div( ) be the divergence operator div acting along the rows of , and t is the transposed tensor of . Besides, the trace of , the sum of the elements on the main diagonal, is denoted by tr( ) and its deviator is defined by d := − 1 tr( )I, where I is the identity matrix of order . Note that the deviator is defined such that tr( d ) = 0. The inner product between the two tensor functions , is given by : = tr( t ), and it follows from the definition that d : = d : d = : d . We consider domains Ω ⊂ R for ∈ {2, 3} with polyhedral boundary Γ := Ω and outward unit normal vector . In terms of notation, for any generic scalar functional space , we denote its vectorial extension by M and its tensor counterpart by M. Furthermore, we employ the standard notation for the space of square-integrable measurable scalar functions L 2 (Ω), along with its norm ‖·‖ 0;Ω , and the Sobolev spaces (Ω) for an integer . We also introduce the spaces: with their norms ‖·‖ div;Ω and ‖·‖ div;Ω , respectively. The norm ‖·‖, with no subscripts, shall be used for any element or operator when it does not cause confusion about the space to which it refers.

The model problem
In this section we present a concise derivation of a model of sedimentation based on the solid phase velocity field. Starting from the equations of conservation of mass and momentum balance of each phase, the derivation of the model equations is mainly based on the ideas of [46] combined with some of the constitutive assumptions made in [16].
We begin by assuming that the sedimentation process occurs on a fixed domain Ω ⊂ R for ∈ {2, 3}, away from sources. Let the volume fraction of solid particles := ( , ) be a function of the position ∈ Ω and time ∈ (0, ], where ∈ R + is the end time. The function is a scalar quantity that takes values between zero and max , the maximum volume fraction. The equations for the conservation of mass of both phases, solid and fluid, are given by where := ( , ) and f := f ( , ), both vector quantities, are the solid and fluid velocity fields, respectively. For the momentum balance of each phase we consider the gravity, viscous and interaction forces, and pressure gradients. Then, neglecting the material derivatives of each phase, the momentum equations are determined by where s := s ( , ) and f := f ( , ) are the pressures of the solid and fluid phase, respectively, s and f are the solid and fluid viscous stress tensors, respectively, is the acceleration of gravity and := (0, 1) t (when = 2) is the upward unit vector. The term present in both equations is the solid-fluid interaction force per unit volume, to be defined later. Adding the momentum equations (2.3) and (2.4) we get the momentum equation of the mixture: where mix := s + f is the viscous stress tensor of the mixture and := ( ) = s + f (1 − ). Now, in order to reduce the number of unknowns we define the excess pore pressure , the hydrostatic pressure h and the effective solid stress function e := e ( ), which is a given function that satisfies where c ∈ (0, max ] is the critical concentration, see [16]. Having these new variables, the pressures of both phases are then given by Moreover, computing the gradient ∇ h = − f , we have The solid-fluid interaction force after neglecting the material derivative of its dynamic part is modeled by a function depending on the relative velocity r := − f , , and h , that is where := ( ) is the Darcy function, which depends on the drag related to the permeability of the material. The solid-fluid interaction force (2.8) is the result of combining both approaches [46] (dynamic part) and [16] (hydrostatic part). Now, in order to obtain a constitutive relation for the relative velocity r , we follow the same idea as in [46]. Replacing and (2.6) into equation (2.4) and assuming that the mixture viscosity of the suspension is orders of magnitude larger than the viscosity of the fluid [46] (which is the case of dense mixtures in mining and wastewater treatment), thus neglecting the tensor f , one arrives at the constitutive relation for r : Defining the volume-average velocity by we obtain div( ) = 0, from which, replacing f = − r = − ( )∇ , it follows that div( − ( )∇ ) = 0, (2.9) where ( ) := (1− ) ( ). We stress here that (2.9) expresses the fact that the divergence of the velocity field is non-zero and varies depending on and , unlike models formulated in terms of the volume average velocity [16,20]. Finally, since f was neglected, we have that mix = s , and we assume that s := s ( , ) is given by s := ( ) ( ), where ( ) := 1 2 (∇ + ∇ t ), and := ( ) is the viscosity function of the mixture which is assumed to be a positive bounded function, i.e., there exist min , max ∈ R + , with min < max , such that Having all the ingredients described here, the main system of equations is where Ω := Ω × (0, ], g( ) := − d − ∇ e ( ) and d = s − f . The first equation is the conservation of mass of the solid phase (2.1), equation (2.11b) is the momentum equation of the mixture (2.5) after replacing the variables and constitutive relations described before, and the third equation corresponds to (2.9). The unknowns of the system are the volume fraction of solid particles , the solid phase velocity and the excess pore pressure . This PDE system is complemented with suitable initial and boundary conditions.

The continuous formulation
In this section we derive a mixed variational formulation for the system (2.11b) and (2.11c) solved for the pair and , at a specific time point ∈ (0, ] for a given function (·, ) ∈ L 2 (Ω), and then apply a fixed-point approach along with the Babuška-Brezzi theory and Banach theorem to address its solvability. For ease of notation and without loss of generality, throughout this section and the subsequent one, the time dependency will be omitted, i.e., the functions (·, ), (·, ) and (·, ) will be simply denoted by , and , respectively, and the same will apply to other variables.

Derivation of the mixed formulation
We begin by introducing some further spaces and related results to be employed in what follows. In particular, we will make use of the decompositions: Note that in light of the second decomposition in (3.1), any element ∈ H(div; Ω) can be uniquely decomposed as = 0 + I, with ∈ H tr (div; Ω) and ∈ R computed by where |Ω| is the measure of the domain Ω. We also define the following spaces: Given a function ∈ L 2 (Ω), the decoupled model equations (2.11b) and (2.11c) are combined with no-slip boundary conditions for and Neumann boundary conditions for , thus yielding We are now going to treat (3.3a) and (3.3b) separately as PDE systems themselves, introducing extra new variables, but coupled through their source terms. Starting from equation (3.3a), we introduce the pseudostress tensor̂︀ := ( ) ( ) − I in Ω, which we seek to belong to H(div; Ω). Considering the space decomposition Then, instead of searching for the pseudostress tensor we focus on the new variable =̂︀ − 0 I which, according to (3.3a), satisfies div( ) = div(̂︀) = −g( ). Note that tr( ) = tr( ) = 0, whence ∈ L 2 tr (Ω). Then, equation (3.4) plus the relations between the new variables , and yield the following system of equations: wherẽ︀ := + 0 is the new unknown pressure. The above system is complemented with an incompressibility hypothesis, which serves as the uniqueness condition for̃︀, namely System (3.5) needs to be solved together with the respective boundary condition from (3.3c), that is = 0 on Γ, and condition (3.6). Note that the pressurẽ︀ can be computed in terms of the new variables and . Indeed, taking trace on both sides of (3.5a), we get which yields̃︀ Furthermore, combining the obtained expression for̃︀ and (3.6) we obtain which, along with the fact that ∈ H tr (div; Ω), implies To derive the weak formulation of (3.5), we test (3.5a) with a tensor ∈ L 2 tr (Ω), which gives ∫︁ Observe that since = d , there holds : = : d = d : . On the other hand, equation (3.5b) can be multiplied by a test function ∈ H tr (div; Ω) and integrated over Ω to obtain where ⟨·, ·⟩ Γ denotes the duality pairing between H −1/2 (Γ) and H 1/2 (Γ). For the last term on the right-hand side of the above equation, we notice that from (3.8) tr( ( )) = (tr( ) +̃︀)/( ( )), so that we get ∫︁ Then, using the boundary condition = 0 on Γ, we arrive at the weak form of (3.5b): The third equation (3.5c) is simply tested against ∈ L 2 (Ω), that is In addition, because of the antisymmetry of , this function is sought in L 2 skew (Ω) and the symmetry of is weakly enforced by ∫︁ Finally, gathering (3.12) and (3.11), we can write Equations (3.9), (3.10) and (3.13) constitute the mixed variational formulation of system (3.5), and if we assume that̃︀ is known, the unknowns of this system are , , and , which will be searched by pairs on the spaces: In this way, giveñ︀ ∈ L 2 0 (Ω) and ∈ L 2 (Ω), the first problem to be solved is: where the bilinear forms and functionals On the other hand, for the weak formulation of (3.3b) we proceed in a similar fashion as for (3.3a) by defining an extra variablẽ︀ := ( )∇̃︀ − , with which we have the system of equations︀ This system is complemented with the boundary conditioñ︀ · = ( ( )∇̃︀ − ) · = 0 on Γ which follows from the second equation in (3.3c). Dividing (3.16a) by ( ), multiplying by a test functioñ︀ ∈ H(div; Ω), and integrating we get ∫︁ Because of the above boundary condition for̃︀, we imposẽ︀ · = 0 on Γ so that the second term on the right-hand side vanishes, as well as the space of test functions becomes Hence, given ∈ L 2 (Ω) and ∈ L 2 (Ω), the second problem reduces to: Find (̃︀,̃︀) ∈ M × such that︀

The fixed-point strategy
Observe that the system (3.14) needs a given pressurẽ︀ ∈ L 2 0 (Ω), which can be obtained from (3.17), while the system (3.17) needs a given velocity ∈ L 2 (Ω), which in turn can be found by solving system (3.14). This coupled nature between the two systems of equations suggests the fixed-point approach to be explained next. In fact, we first let S : → L 2 (Ω) be the operator defined for each fixed function ∈ L 2 (Ω) by where (( , ), ( , )) ∈ X × V is the solution of (3.14) for the giveñ︀ ∈ and ∈ L 2 (Ω). In turn, we let︀ S : L 2 (Ω) → be the operator defined for each fixed function ∈ L 2 (Ω) bỹ︀ where (̃︀,̃︀) ∈ M × is the solution of (3.17) for the given ∈ L 2 (Ω) and ∈ L 2 (Ω). The well-posedness of S and̃︀ S , or equivalently the uniqueness of solutions of both systems (3.14) and (3.17), respectively, is shown next in Section 3.3. Now, defining the operator T =̃︀ S ∘ S : → , we readily realize that solving the coupled problem (3.14)-(3.17) is equivalent to finding a fixed-point of T , that is: find̃︀ ∈ such that T (̃︀) =̃︀. (3.19) We end this section by remarking that the operators S and̃︀ S are affine and linear, respectively, which follows from the properties of ,̃︀ and̃︀ , , and hence T becomes affine as well.
Finally, using from the definitions of ,̃︀ and (cf. (3.15)) that we arrive at (3.21) and finish the proof.
From now on for the second system (3.17), we assume that max < 1 and that there exist min , max ∈ R + such that for < max min ≤ ( ) ≤ max . (3.23) The following lemma employs the classical Babuška-Brezzi theory to show the existence and uniqueness of solutions for (3.17) in the space M × . Proof. We begin by observing that both bilinear forms̃︀ a and̃︀ b, and̃︀ , are bounded, namely: which shows that̃︀ a is elliptic in kernel(̃︀ b). On the other hand, the bilinear form̃︀ satisfies the following inf-sup condition: with 3 > 0, see [42] for further details. Then, a straightforward application of Theorem 2.3 from [42], and the fact that ‖̃︀ , ‖ M ′ ≤ ‖( ( )) −1 ‖ 0;Ω (cf. (3.24)) conclude the proof.
Note that Lemmas 3.1 and 3.2 imply the well-posedness of the operators S and̃︀ S , respectively, and therefore their composition, the operator T , is well defined.

Well-posedness of the fixed-point equation
The aim of this section is to apply the Banach fixed-point theorem to show that (3.19) has a unique solution. To this end, for a given > 0, we define the ball B := { ∈ : ‖ ‖ 0;Ω ≤ }. The next lemma shows that the operator T maps B into itself for sufficiently small data. Lemma 3.3. Given > 0 and ∈ L 2 (Ω) such that < max on Ω, we assumẽ︀

25)
where and̃︀ are the constants determined by Lemmas 3.1 and 3.2, respectively. Then, the inclusion T (B ) ⊆ B holds.
Continuing, we now show that both operators S and̃︀ S are Lipschitz continuous.
Lemma 3.4. Let ∈ L 2 (Ω) such that < max on Ω. Then, the following inequalities hold Proof. To prove the first inequality, we make use that S is affine. In this way, giveñ︀ 1 ,̃︀ 2 ∈ L 2 0 (Ω), S (̃︀ 1 ) − S (̃︀ 2 ) returns the solution of system (3.14) whose right hand side of (3.14b) is evaluated at̃︀ :=̃︀ 1 −̃︀ 2 , and the right hand side of (3.14c) is zero. Then, replacing = 0 in Lemma 3.1, we get Analogously, using now the linearity of̃︀ S and Lemma 3.2 instead, we find that We are now in condition to introduce the theorem that states the existence and uniqueness of weak solutions of the coupled system (3.14)-(3.17), which is the weak version of (3.3).
Theorem 3.5. Let ∈ L 2 (Ω) such that < max in Ω, and let > 0. In addition, let ,̃︀ be the constants determined by Lemmas 3.1 and 3.2, respectively, and assume that the data satisfy (3.25) and︀ min min < 1. (3.26) Then, the coupled system (3.14)-(3.17) has a unique solution (( , ), ( , )) ∈ X × V and (̃︀,̃︀) ∈ M × , with︀ ∈ B , and there hold Proof. From the assumption (3.25), Lemma 3.3 ensures that T (B ) ⊆ B . Then, making use of Lemma 3.4 and the properties of̃︀ S , S , and T , for each̃︀ 1 ,̃︀ 2 ∈ we have and since min min >̃︀ , we get that T is a contraction. The uniqueness of solution for the equation T (̃︀) =̃︀, or equivalently of the coupled system (3.14)-(3.17) with̃︀ ∈ B , is a result of the application of the Banach fixed-point theorem. The a priori estimates follow from Lemmas 3.1 and 3.2.
At this point we remark that (3.26) constitutes a technical restriction of the analysis, which, according to the lower bounds from (2.10) and (3.23), requires, on one hand, that both, the viscosity function of the mixture and the Darcy function-dependent coefficient , attain sufficiently large values. In turn, the boundedness constants (cf. Lem. 3.1) and̃︀ (cf. Lem. 3.2) need to be small enough in order to guarantee (3.26). In this way, and additionally to being certainly a somehow restrictive assumption from the physical point of view, mainly because of the above constraints on and , we notice that (3.26) is actually not verifiable in practice since the aforementioned constants are usually not known explicitly.

Discrete formulation
In this section we introduce the Galerkin scheme associated with the coupled system given by (3.14) and (3.17), and proceed analogously to the analysis developed in Section 3 to derive its well-posedness.

Preliminaries and setting of the scheme
We begin by letting { ℎ } ℎ>0 be a regular family of triangulations of Ω made up of triangles (when = 2) or tetrahedra (when = 3) of diameter ℎ , where ℎ stands for both a sub-index of each triangulation ℎ , and for the largest diameter, that is ℎ := max{ℎ : ∈ ℎ }. Next, given an integer ℓ ≥ 0 and a subset Λ of R , we denote by ℓ (Λ) the space of polynomials of total degree less than or equal to ℓ defined on , ℓ (Λ) := [ ℓ (Λ)] and P ℓ (Λ) := [ ℓ (Λ)] × . We also define¯ℓ the set of polynomials of degree equal to ℓ. Then, denoting by a generic vector of R , the local Raviart-Thomas space of order ℓ is defined for each ∈ ℎ by RT ℓ ( ) := ℓ ( ) ⊕¯ℓ( ) .
In turn, the local bubble space is defined by where is the bubble function on , which is defined as the product of its + 1 barycentric coordinates, curl( ) = (︁ , − )︁ for = 2 and : → R, and curl( ) = ∇ × for = 3 and : → R 3 . Furthermore, we set the global spaces: and their corresponding tensor versions  Hence the discrete subspaces for X × V are X ℎ := H ℎ × H ℎ and V ℎ := H ℎ × H ℎ . We observe that M ℎ contains the multiples of I, and in consequence we can equivalently define H ℎ by Moreover, in light of the fact that div(M ℎ ) ⊆ H ℎ , the discrete kernel of the operator ℬ 2 reduces to 0,ℎ := The subspaces for the approximation of (3.17) are chosen by Hence, giveñ︀ ℎ ∈ ℎ and ∈ L 2 (Ω), the discrete formulation of (3.14) reads: find (( ℎ , ℎ ), ( ℎ , ℎ )) ∈ X ℎ ×V ℎ such that In turn, given ℎ ∈ H ℎ and ∈ L 2 (Ω), the discrete formulation of system (3.17) is given by: To state the discrete version of (3.19) we first define the discrete counterparts of the operators S ,̃︀ S , and T , for which we consider here that ∈ L 2 (Ω) is a given function. Let S ,ℎ : ℎ → H ℎ be the operator defined by where (( ℎ , ℎ ), ( ℎ , ℎ )) ∈ X ℎ × V ℎ is the unique solution of (4.3) for the giveñ︀ ℎ ∈ ℎ and ∈ L 2 (Ω). Similarly, we let̃︀ S ,ℎ : H ℎ → ℎ be the operator defined bỹ︀ where (̃︀ ℎ ,̃︀ ℎ ) ∈ M ℎ × ℎ is the unique solution of (4.4) for the given ℎ ∈ H ℎ and ∈ L 2 (Ω). The wellposedness of S ,ℎ and̃︀ S ,ℎ , which is equivalent to proving the uniqueness of solutions of (4.3) and (4.4), respectively, is shown in Section 4.2. Finally, we define the operator T ,ℎ :=̃︀ S ,ℎ ∘ S ,ℎ : → , and realize that the coupled problem (4.3) and (4.4) is equivalent to finding a fixed-point of T ,ℎ , that is̃︀ ℎ ∈ such that Analogously to the continuous case, S ,ℎ ,̃︀ S ,ℎ , and T ,ℎ are affine, linear, and affine, respectively.

Solvability of the discrete problem
The solvability analyses of the discrete problems (4.3) and (4.4), and of the fixed-point equation (4.5), are carried out using arguments analogous to those used for the continuous case. In particular, the next two results are the discrete versions of Lemmas 3.1 and 3.2, respectively. Lemma 4.1. Giveñ︀ ℎ ∈ ℎ and ∈ L 2 (Ω), there exists a unique (( ℎ , ℎ ), ( ℎ , ℎ )) ∈ X ℎ × V ℎ solution of (4.3), and hence one can define S ,ℎ (̃︀ ℎ ) := ℎ . Moreover, there exists a positive constant d > 0, depending on the dimension , min , and max , and thus independent of ℎ, such that Proof. The Lipschitz continuity and strong monotonicity of were shown in the proof of Lemma 3.1. On the other hand, since [ ( ℎ ), ℎ ] ≥ 0 for all ℎ ∈ H ℎ , is clearly positive semi-definite in ℎ,0 , the discrete kernel of ℬ 2 . In addition, given ℎ ∈ 0,ℎ , we have that div( ℎ ) = 0 and ‖ d ℎ ‖ 0;Ω ≥ 0 ‖ ℎ ‖ div;Ω ( [42], Lem. 2.3), and since d ℎ ∈ H ℎ , we find that which shows that the operator ℬ 1 satisfies the discrete inf-sup condition. In turn, the same condition for ℬ 2 has been shown in previous works, e.g., [44]. Consequently, since the scheme satisfies the hypotheses of Theorem 3.4 from [43], the proof is concluded.
We notice here that essentially the same remark stated at the end of Section 3.4 regarding (3.26), is valid for (4.6). Nevertheless, we stress in advance that the numerical results to be reported later on in Section 6 are not affected at all by the eventual lack of verification of this assumption, thus suggesting that the latter is just a limitation of the analysis rather than the method.

A priori error analysis
be the bilinear forms arising from adding the left hand sides of (3.14) and (3.17), respectively. Then, introducing the generic notation the pairs of continuous and discrete schemes (3.14)-(4.3) and (3.17)-(4.4) can be rewritten as respectively. Thus, as a consequence of Lemmas 4.1 and 4.2, we deduce the existence of positive constants A and̃︀ A , independent of ℎ, with which A and̃︀ A satisfy global discrete inf-sup conditions, namely and Hereafter, given a subspace ℎ of a generic normed space ( , ‖ · ‖ ), we set for each ∈ dist( , ℎ ) := inf The following lemma states the Céa estimate for our coupled problem.
Lemma 4.6. Let ∈ L 2 (Ω) such that < max in Ω, and assume that the data satisfies √ min A Then, there exist positive constants 1 , 2 , depending on ‖A ‖, A , and ‖̃︀ A ‖,̃︀ A , respectively, and thus independent of ℎ, such that Proof. It proceeds similarly as in [44]. In fact, bearing in mind (4.9) and (4.10), straightforward applications of the classical Strang estimate to the pairs (4.7) and (4.8) yield and respectively. In turn, it follows from (3.22) and (2.10) that whereas the bound for ‖̃︀ , ‖ (cf. (3.24)) and (3.23) give In this way, employing (4.15) and (4.16) back into (4.13) and (4.14), respectively, and then adding the resulting inequalities, and using (4.11), we obtain (4.12) with We point out now that similar remarks to those provided at the end of Sections 3.4 and 4.2 hold here for (4.11). We omit further details and refer to those statements.
We can now state the rates of convergence of the Galerkin scheme (4.3) and (4.4) with the finite element subspaces (4.1) and (4.2).

Fully coupled scheme
In this section we address the approximation of the transport equation (2.11a) for the volume fraction of solid particles , and the coupling to the systems (4.3) and (4.4). Since (2.11a) is of hyperbolic type, we propose a vertex-centered finite volume scheme based on an upwind numerical flux for its approximation. In what follows, we let ℎ be the set of nodes of the triangulation ℎ (primal mesh). We introduce the so-called dual mesh ⋆ ℎ which consists of the non-overlapping control volumes ⋆ surrounding the node ∈ ℎ . The control volumes on the dual mesh ⋆ ℎ are built by connecting the barycenters of each triangle ∈ ℎ with the middle point of each edge of , see Figure 1. Note that each node in the primal mesh is related to one control volume in the dual mesh. We denote by ℰ ℎ ( ⋆ ) the set of inner edges of ⋆ ∈ ⋆ ℎ , so that ∈ ℰ ℎ ( ⋆ ) if ⊂ ⋆ and there exists ∈ ℎ such that is the piece of segment connecting and the middle point of an edge of (see the blue dash-dotted segments in Fig. 1). In turn, we define the space of piecewise constant functions over the dual mesh To introduce the finite volume method, we consider the integral form of equation (2.11a) over each control volume in the dual mesh ⋆ ℎ , this is ∫︁ Then, the sought unknown is approximated by ℎ ( ) ∈ ℎ , i.e., for each ∈ (0, ], is a piecewise constant function on ⋆ ℎ . Furthermore, given a time step ∆ > 0 and ∈ N, we make use of the superscript to denote the evaluation of each function at the time = ∆ , e.g., ℎ = ℎ ( ), ℎ = ℎ ( ), and the analogous notation is used for other functions. Next, approximating the time derivative in (5.1) by forward Euler and the numerical flux by an upwind approach we obtain where ℎ, ⋆ = ℎ | ⋆ , is the normal to the edge ∈ ℰ ℎ ( ⋆ ), ⋆ ∈ ⋆ ℎ denotes the adjacent control volume such that ⊆¯⋆ ∩¯⋆ and ⋆ ̸ = ⋆ , and the upwind operator is defined by Upw( ; , ) := max{ , 0} + min{ , 0} , , , ∈ R.
The velocity function ℎ is determined by solving (4.4) and (4.3) (the fixed-point Eq. (4.5)) at the time with given by the piecewise constant function that for each ∈ ℎ is computed by Note the difference between ℎ and , the first one belongs to the subspace ℎ (defined over the dual mesh), while the second one is in 0 (Ω) (defined over the primal mesh). For the case of a piecewise linear (or constant) approximation ℎ , the integral over the edge in (5.2) can be computed with the Gauss quadrature rule to obtain where | | and | ⋆ | are the measures of the edge and control volume ⋆ , respectively, and¯is the middle point of . The discrete approximation of (5.3) needs to be solved for a sufficiently small time step ∆ .
where coeff ( ) is the vector of coefficients of̃︀ ℎ at the iteration ∈ N, | · | is the Euclidean norm of R dim( ℎ ) , and tol is a positive number. For the rest of the examples, equation (4.5) is solved by the Newton-Raphson algorithm with null initial guess, with absolute and relative tolerance of 10 −10 , and the solution of tangent systems resulting from the linearization is carried out with the multifrontal massively parallel sparse direct solver MUMPS [7]. In addition, the integral condition (3.6) is added to the system

Example 1. Accuracy verification
We test the order of accuracy of the numerical scheme used to approximate the decoupled system (4.3) and (4.4) in the unit square Ω = (0, 1) 2 for a given smooth volume fraction < max . For this purpose, we set 0 = 1 and 0 = 1, and consider Furthermore, we use the manufactured exact solution: for which the right hand sides of (4.3c) and (4.4b) have to be modified by adding the terms respectively. The errors between each component of the numerical solution ℎ , ℎ , ℎ , ℎ ,̃︀ ℎ ,̃︀ ℎ , computed with a meshsize ℎ, and the exact manufactured solution, are given by The degrees of freedom for the system (4.3) and (4.4), denoted by dof1 and dof2, respectively, and the errors and rates of convergence produced by successive mesh refinements, are reported in Table 1 (for ℓ = 0) and 2 (for ℓ = 1). The number of iterations made by the fixed-point algorithm with a tolerance of tol = 10 −6 was 6 for all the examples, for ℓ = 0 and ℓ = 1. From Table 1, we observe that for all variables the error decreases as the mesh size is smaller and the rates of convergence tend to 1, with the exception of whose rate of convergence is greater than 1. For the case of polynomial degrees ℓ = 1 in Table 2, the errors behave as expected with rates of convergence approaching 2 as Theorem 4.7 states. Hence, both examples show that for ℓ = 0 and ℓ = 1 the errors produced with the scheme (4.3) and (4.4), and the fixed-point strategy, are of order (ℎ ℓ+1 ), the optimum established in Theorem 4.7. Additionally, in Figure 2 we show the numerical solution of the six unknowns produced in the accuracy test. Note that the order of accuracy shown in Tables 1 and 2 does not take into account the approximation of the transport equation.

Examples 2 and 3
For the second and third example, we simulate the transient case of two-dimensional batch sedimentation of flocculated particles in a channel with inclined walls. The domain considered in both simulations is a truncated cone of height 2 m and 6 m wide, and represents the vertical cross-section of a channel. In Example 2 we start the simulation with an initial concentration where below the level 2 = 1.5 m the concentration is homogeneous and above it there is only clear water. For Example 3 we consider an initial condition whose concentration of solid particles is located above the line otherwise.
In Figure 3, left column, we show the simulation of Example 2 at four time points = 4, 6, 8 and 10 s. We observe that as soon as time evolves, the sedimentation process takes place, and the suspended solid particles move downward concentrating at the bottom of the channel, as expected. In accordance with the no-slip boundary condition, the solid-fluid interface moves faster at the middle of the channel and the concentration of solid particles increases towards the side walls and bottom. Furthermore, on each time iteration, the conservation of mass is fulfilled. In addition, in the second column of Figure 3, we show the simulation of Example 3 at the time points = 5, 15, 20 and 35 s. We observe the downward falling of the solid particles which tend to concentrate at the boundaries and move towards the bottom. A layer of concentrated particles is clearly divided from a region of low concentration that shows the path that particles have followed in the sedimentation process. The movement of the solid particles is not only vertical but also directed to the right due to accumulation of particles at the left side wall. In this simulation the conservation of mass is also satisfied.

Examples 4 and 5
In the fourth and fifth example, we perform the batch sedimentation on an inclined channel of rectangular cross-section with an angle of inclination of 30 ∘ degrees with respect to the vertical. The length and width of the rectangular domain are 10 and 2 m, respectively. In Example 4 the initial condition chosen is and for Example 5 we consider an initial concentration varying with the width of the rectangle given by In Figure 3, first row, we show the simulation of Example 4 at the time points = 6, 9 and 12 s. We observe that the solid particles tend to settle towards the down side wall (right wall) and bottom as expected. However, at the  top side wall (left wall) near the line 2 = 8 the solid particles do not separate from the clear water immediately, but begin to settle at a slower pace, which can be explained as a consequence of the no-slip boundary condition. Example 5 is shown in the bottom row of Figure 3, for the time points = 6, 15 and 20 s. In this case the initial concentration is zero at the left side wall and increases continuously up to 0.1 at the right side wall. As in the previous example the concentration increases towards the right side and bottom wall but due to the

Discussion
We have introduced a model of multidimensional sedimentation based on the solid phase velocity which can be seen as a simplified version of [46]. From the model consisting of three unknowns and equations, we have studied the solutions of the decoupled velocity-pressure equations by using mixed finite element approaches. The uniqueness of the decoupled system is shown based on the assumption that the function ( ) = (1 − ) ( ) is bounded and away from zero. Despite the fact that the analysis of existence and uniqueness of solutions relies on this condition, the case when ( ) = 0 reduces to the classical Stokes system which has been widely analyzed in the literature [42]. One of the main advantages of having a mixed formulation is that extra variables, such as the pseudostress and vorticity tensors, can be recovered as a component of the solution of the problem. Such variables may have practical interest in the study of, for instance, the viscosity of mixtures [50]. Although the finite element subspaces for system (4.3) are based on PEERS finite elements [9], another suitable election could be the Arnold-Falk-Winther subspaces [10], which are also stable.
A finite-volume method making use of an upwind flux approximation over a dual mesh is proposed to get conservation of mass for the solid volume fraction when coupling to the mixed finite element scheme for the velocity, pressure and additional variables. From the numerical examples we conclude that the coupled numerical scheme performs well, it is stable under sufficiently small time step and gives physically relevant numerical solutions (cf., e.g., [14,20]). Despite the simple geometries considered in the five examples, the scheme can be implemented for domains with more involved geometries and even in three-dimensions. Another advantage of our model, as compared with, for example, the classical -model [16], is that the transport equation is linear in and a simple upwind scheme is sufficient to obtain physically relevant approximate solutions.
The model equations, analysis and numerical scheme can be extended to include, for example, a more general viscous stress tensor (eventually non-linear in ), a non-zero boundary condition for the velocity field or to consider constitutive functions and being able to assume zero values. Further studies can be done to compare the model presented here with related systems, and also with problems such as the one in [16].