FINITE ELEMENT METHODS FOR THE DARCY-FORCHHEIMER PROBLEM COUPLED WITH THE CONVECTION-DIFFUSION-REACTION PROBLEM

In this article, we consider the convection-diffusion-reaction problem coupled the DarcyForchheimer problem by a non-linear external force depending on the concentration. We establish existence of a solution by using a Galerkin method and we prove uniqueness. We introduce and analyse a numerical scheme based on the finite element method. An optimal a priori error estimate is then derived for each numerical scheme. Numerical investigation are performed to confirm the theoretical accuracy of the discretization.


Introduction
This work studies the convection-diffusion-reaction equation coupled with Darcy-Forchheimer problem. The system of equations is ( ) where Ω ⊂ IR , = 2, 3, is a bounded simply-connected open domain, having a Lipschitz-continuous boundary Γ with an outer unit normal n. The unknowns are the velocity u, the pressure and the concentration of the fluid. |.| denotes the Euclidean norm, |u| 2 = u · u. The parameters , and represent the density of the fluid, its viscosity and its dynamic viscosity, respectively. is also referred as Forchheimer number when it is a scalar positive constant. The diffusion coefficient and the parameter 0 are positive constants. The function Keywords and phrases. Darcy-Forchheimer problem, convection-diffusion-reaction equation, finite element method, a priori error estimates.
f represents an external force that depends on the concentration and the function represents an external concentration source. is the permeability tensor, assumed to be uniformly positive definite and bounded such that there exist two positive real numbers and such that It should be noted that should be smaller than the smallest eigenvalue of −1 over Ω and could be very large.
To simplify, a homogeneous Dirichlet boundary condition is prescribed on the concentration , but the present analysis can be easily extended to a non-homogeneous boundary condition.
System ( ) couples the Darcy-Forchheimer system with the convection-diffusion-reaction equation satisfied by the the concentration of the fluid. The same system can also couple the Darcy-Forchheimer system with the heat equation satisfied by the temperature of the fluid, it suffices to set 0 = 0 and replace by .
Darcy's law (see [31] and [37] for instance for the theoretical derivation) describes the creeping flow of Newtonian fluids in porous media. It is simply the first equation of system ( ) without the non-linear term |u|u and where the function f may depend on the concentration of the fluid. Forchheimer [19] showed experimentally that when the velocity is higher and the porosity is nonuniform, Darcy's law becomes inadequate. He proposed the Darcy-Forchheimer equation which is the first equation of the system ( ). A theoretical derivation of Forchheimer's law can be found in [29]. Multiple theoretical and numerical studies of the Darcy-Forchheimer system were performed and among others we mention [24,26,28,32,33]. Many numerical investigations are performed and show the importance of the Darcy-Forchheimer equation compared with Darcy equation (see for instance [32] and the references inside).
For the coupled problem of Darcy's law with the heat equation, we can refer to [6] where the coupled problem is treated using the spectral method. The authors in [4] and [14] considered the same stationary system but coupled with a nonlinear viscosity that depends on the temperature. In [15], the authors derived an optimal a posteriori error estimate for each of the numerical schemes proposed in [4]. We can also refer to [3] where the authors used a vertex-centred finite volume method to discretize the coupled system. Furthermore, for the time-dependent convection-diffusion-reaction equation coupled with Darcy's equation, we refer to [9,10] where the authors established the corresponding a priori and a posteriori errors.
The coupling system ( ) has many physical applications for the Darcy-Forchheimer mixed convection case [36]. In this case, the Darcy-Forchheimer system is coupled with the concentration of the fluid with an external force f . We mention that the work of [6] use the spectral method to treat a coupled system very close to ( ) where the Darcy-Forchheimer equation is replaced by the Darcy's one. It turns out that the non-linear term appearing in the first equation of ( ) makes the treatment of the coupled system more complex.
We first derive an equivalent variational formulation to ( ) and we show the existence of a solution. The uniqueness can be reached under additional constraint on the concentration (see Condition (2.27)). Then, we discretize the system by using the finite element method and we show the existence and uniqueness of the corresponding solution. Later, we establish the a priori error estimate between the exact and numerical solutions under the condition of smallness of the concentration in the fluid. In order to compute the solution, we introduce an iterative scheme and we study the corresponding convergence. Finally, numerical investigations are performed to validate the theoretical results.
The outline of the paper is as follows: -Section 2 is devoted to the continuous problem and the analysis of the corresponding variational formulation.
-In Section 3, we introduce the discrete problems, recall their main properties, study their a priori errors and derive optimal estimates. -In Section 4, we introduce an iterative algorithm and prove its convergence.
-Numerical results validating the convergence analysis are presented in Section 5.

Notation
Let (Ω) be the space of functions having a compact support in Ω with continuous derivatives of all orders in Ω. Let = ( 1 , 2 , . . . ) be a -uplet of non-negative integers, set | | = ∑︀

=1
, and define the partial derivative by Then, for any positive integer and number ≥ 1, we recall the classical Sobolev space [2,30] , When = 2, this space is the Hilbert space (Ω). The definitions of these spaces are extended straightforwardly to vectors, with the same notation, but with the following modification for the norms in the non-Hilbert case. Let v be a vector valued function; we set We shall often use the following Sobolev embeddings: for any real number ≥ 1 when = 2, or 1 ≤ ≤ 2 −2 when ≥ 3, there exist constants and 0 such that

Variational formulation
In this section, we introduce the variational formulation corresponding to the problem ( ). We assume that the volumic and boundary sources verify the following conditions: We assume that f and verify: (1) f can be written as follows: (Ω) and f 1 is Lipschitz-continuous with constant f1 and satisfies the inequality where f1 is a strictly positive constant. (2) ∈ 2 (Ω).

Remark 2.2.
For the physical interpretation of the choice of the external force f , we refer to [6,7,13,20]. In fact, in [6] page 3, they considered and justified the choice of f ( , ) = f 1 ( ) (with f 0 = 0) with the corresponding properties. For the generalization, we added the function f 0 to get (2.1).
It follows from the nonlinear term in the system ( ) that the velocity u and the test function v must belong to 3 (Ω) ; then, the gradient of the pressure must belong to 3 2 (Ω) . Furthermore, the concentration must be in 1 0 (Ω). Thus, we introduce the spaces Furthermore, we recall the following inf-sup condition between and [24], With these assumptions on the sources, we introduce the following variational formulation associated to problem ( ): Equivalence between ( ) and ( ) in the sense of distribution follows readily from the validity of Green's formula: and the fact that for further details, we refer to [24].
Thus, problem ( ) can be rewritten as a function of the single unknown . Indeed, for a given , let (u( ), ( )) be the solution of problem (2.12). Then, problem ( ) is equivalent to the following reduced formulation: find ∈ such that Before proving that problem (2.15) admits a solution, we will show the following intermediate lemma: Lemma 2.5. Under Assumption 2.1, let ( ) ≥1 be a sequence of functions in 2 (Ω) that converges strongly to in 2 (Ω). Then, the sequence (u( )) ≥1 converges strongly to u( ) in and the sequence ( ( )) ≥1 converges weakly to ( ) in .
Finally, we have to treat the convergence of the pressure. Since u( ) is a solution of problem (2.12), we use the inf-sup condition (2.10) to deduce the existence of ( ) such that (u( ), ( )) is the solution of problem (2.12).
We deduce from problem (2.12) that for all v ∈ The strong convergence of u( ) to u( ) in 3 (Ω) , of f (., ) to f (., ) in  The next theorem shows the existence of at least one solution to the problem ( ). Theorem 2.6. Under Assumption 2.1, problem ( ) admits a solution in × × . Furthermore, each solution (u, , ) of ( ) satisfies the following bounds: (2.20) Proof. We propose to construct a solution of ( ) by Galerkin's method. As 1 0 (Ω) is separable, it has a countable basis ( ) ≥1 . Let Θ be the space spanned by the first basis functions, ( ) 1≤ ≤ . Problem (2.13) is discretized in Θ by the square system of nonlinear equations: find = ∑︁ =1 ∈ Θ solution of: where (u( ), ( )) solves (2.12) with = . Now, given ∈ Θ , we introduce the auxiliary problem: find Φ( ) ∈ Θ such that, for all ∈ Θ , we have In fact, for each , the left hand side of (2.22) is a bilinear map of the form (∇., ∇.) and the right hand side is a linear map with respect to . Relation (2.22) defines a continuous mapping from Θ into Θ , due to the fact that Θ is a finite dimension space and Lemma 2.5. By taking = , we get )︁ .
In other words, we have Therefore, Brouwer's Fixed-Point Theorem implies immediately the existence of at least one solution to the problem (2.21).
Let be a solution of problem (2.21), satisfying ∀ ∈ Θ , By taking = in the last equation, we get immediately the bound The last uniform bound implies that, up to a subsequence, ( ) converges weakly to a function in 1 0 (Ω). Therefore, it converges strongly in (Ω), for any < 2 − 2 , and it follows from Lemma 2.5 that (u( ), ( )) converges weakly to (u( ), ( )) in × , and (u( )) converges strongly to u( ) in 3 (Ω) . Now, we freeze the index in (2.21), and let tends to infinity. The weak convergence of ( ) to in 1 0 (Ω), and the strong convergence of (u( )) to u( ) in 3 (Ω) allow us to deduce that is a solution of the following problem: (2.23) From this system and the density of the basis in 1 0 (Ω), we infer that is a solution of problem (2.15). The first bound in (2.20) can be straightly obtained by taking = in the last equation in (2.15). The first and second bounds in (2.20) can be deduced from the inequalities (2.14) and Assumption 2.1.
Theorem 2.7. Assume that Ω is of class 1,1 . Let (u, , ) be a solution of problem ( ). If ∈ ∞ (Ω), then the concentration is in ∞ (Ω) and satisfies the following bound: Proof. Let (u, , ) be a solution of problem ( ), then the velocity u ∈ . Using the fact that is separable and that the space is dense in (see for instance [18], Lem. 10.8), there exists a sequence (u ) ∈N in ∞ which converges strongly to u in 3 (Ω) when tends to +∞. Now, for each ∈ N, let ∈ 1 0 (Ω) the unique solution to the following problem: The elliptic regularity (see [21]) allows us to get that ∈ 2, (Ω) ∩ 1 0 (Ω), for all ≥ 1. Multiplying the first equation of (2.24) by 2 +1 and integrating by parts give, By remarking that the first term of the left hand side of the previous equation is non-negative and by applying Holder's inequality for the right-hand-side with the conjugate exponents = 2 + 2 2 + 1 and = 2 + 2, we get the following inequality: The next step shows that converges strongly to ∈ 1 0 (Ω). In order to prove it, we start by subtracting the third equation of problem ( ) from first equation of (2.24) to get, for all ∈ 1 0 (Ω), By taking = − and using the antisymmetric property, the estimate (2.20) we obtain which gives the strong convergence of to in 1 0 (Ω). As is uniformly bounded in 2 +2 (Ω), we can extract a subsequence still denoted by such that converges weakly in 2 +2 (Ω) to some function ℎ satisfying (2.25). The strong convergence of to in 1 (Ω) and the uniqueness of the limit allows us to deduce that ℎ = in 2 +1 2 +2 (Ω) and we get Thus ℎ = ∈ 2 +2 (Ω) and finally, taking the limit in (2.26) as → ∞, we get the desired result.

27)
then, the solution of the problem ( ) is unique.
Corollary 2.9. Under Assumption 2.1 and Theorems 2.7 and 2.8, if the data satisfies the following smallness condition

Discretization
From now on, we assume that Ω is a polygon when = 2 or polyhedron when = 3, so it can be completely meshed. For the space discretization, we consider a regular (see Ciarlet [11]) family of triangulations ( ℎ ) ℎ of Ω which is a set of closed non degenerate triangles for = 2 or tetrahedra for = 3, called elements, satisfying -for each ℎ,Ω is the union of all elements of ℎ ; -the intersection of two distinct elements of ℎ is either empty, a common vertex, or an entire common edge (or face when = 3); -the ratio of the diameter ℎ of an element ∈ ℎ to the diameter of its inscribed circle when = 2 or ball when = 3 is bounded by a constant independent of ℎ, that is, there exists a strictly positive constant independent of ℎ such that, As usual, ℎ denotes the maximal diameter of all elements of ℎ . To define the finite element functions, let be a non-negative integer. For each in ℎ , we denote by P ( ) the space of restrictions to of polynomials in variables and total degree at most , with a similar notation on the faces or edges of . For every edge (when = 2) or face (when = 3) of the mesh ℎ , we denote by ℎ the diameter of . We shall use the following inverse inequality [15]: for any dimension , there exists a constant such that for any polynomial function ℎ of degree on , The constant depends on the regularity parameter of (3.1), but for the sake of simplicity this is not indicated.
Let ℎ ⊂ , ℎ ⊂ , and ℎ ⊂ be the discrete spaces corresponding to the velocity, the pressure and the concentration. We assume that ℎ and ℎ satisfy the following inf-sup condition: where 2 is a strictly positive constant independent of ℎ. Problem ( ) can be discretized as following: In the following, we will introduce the finite dimension spaces ℎ , ℎ and ℎ . Let be an element of ℎ with vertices , 1 ≤ ≤ + 1, and corresponding barycentric coordinates . We denote by ∈ P +1 ( ) the basic bubble function: We observe that (x) = 0 on and that (x) > 0 in the interior of .
We introduce the following discrete spaces: In this case, for the inf-sup condition (3.3), we refer to [23]. We shall use the following results: (1) For the concentration: there exists an approximation operator (when = 2, see Bernardi and Girault [5] or Clément [12]; when = 2 or = 3, see Scott and Zhang [35]), ℎ in ℒ( 1, (Ω); ℎ ) such that for all in ℎ , = 0, 1, = 0, 1, and all ≥ 1, where ∆ is the macro element containing the values of used in defining ℎ ( ). (2) For the velocity: we introduce a variant of ℎ denoted by ℱ ℎ (see [4] and [22]) which is stable in 3 (Ω) : such that ℱ ℎ (v) ∈ ℎ when div v = 0, and satisfies (3.7). (3) For the pressure: as ℎ contains all constants, an easy modification of ℎ yields an operator ℎ ∈ ℒ( 1, (Ω) ∩ 2 0 (Ω); ℎ ) (see for instance Abboud et al. [1]), satisfying (3.7). Indeed, ℎ can be constructed as follows: Existence of a solution of ( ℎ ) is derived by duplicating the steps of the previous section concerning the existence of a solution of problem ( ). First ( ℎ ) is split as in the previous section, i.e., find ℎ ∈ ℎ such that: (3.10) For each ℎ ∈ ℎ , an easy finite-dimensional variant of the argument of Theorem 2.4 allows one to prove that the scheme (3.10) has a unique solution (u ℎ ( ℎ ), ℎ ( ℎ )) ∈ ℎ × ℎ , and this solution satisfies the a priori estimates similar to (2.14): . (3.11) We address now the existence of at least one solution of the problem (3.9) written with the only variable ℎ . For this purpose, we apply Brouwer's Fixed-Point Theorem. Indeed, we introduce the following map: for a given ℎ ∈ ℎ , find Φ( ℎ ) ∈ ℎ such that: This last relation defines a mapping from ℎ into itself, and we easily derive its continuity. By taking ℎ = ℎ , we get In other words, we have (∇Φ( ℎ ), ∇ ℎ ) ≥ 0, for all ℎ ∈ ℎ such that The Brouwer's Fixed-Point Theorem implies immediately the existence of at least one solution of the problem (3.9). Hence, problem ( ℎ ) admits at least one solution (u ℎ , ℎ , ℎ ) ∈ ℎ × ℎ × ℎ . Furthermore, by taking ℎ = ℎ in the last equation of ( ℎ ) gives, in addition to inequality (3.11), the following bound: (3.12) Finally, uniqueness follows easily since ℎ belongs to ∞ (Ω). This is summed up in the following existence and uniqueness theorems.
The next step consists to show that u ℎ converges strongly toū in 3 (Ω) . In order to prove it, we start by taking v = v ℎ in (2.12) and subtracting from (3.10), we have (3.26) By inserting (ℱ ℎ (ū)) and taking v ℎ = u ℎ − ℱ ℎ (ū), we get (3.27) The monotonicity of allows us to obtain, (3.28) We pass to the limit in the previous equation. We deduce from the strong convergence of (ℱ ℎ (ū)) to (ū) and of f (., ℎ ) to f (.,¯) that the first and last terms of the right-hand-side of the previous inequality tend to 0. Furthermore, the weak convergence of u ℎ toū, and the strong convergence of ℱ ℎ (ū) toū imply the convergence of the second term of the right hand side of the previous inequality to 0. Thus, u ℎ converges strongly toū in 3 (Ω) .
To finish the proof, it remains to show thath = ∇¯, which can be easily obtained by passing to the limit in (3.26), and by using the strong convergence of u ℎ toū in 3 (Ω) , and the uniqueness of the weak limit. Proof. We have proved in Proposition 3.4 that (ū,¯,¯) solves the first two equations of problem ( ). It remains to show that (ū,¯) solves the third equation of problem ( ). We consider the third equation of problem ( ℎ ). By taking ℎ = ℎ for a regular ∈ (Ω) (taking into account the density of (Ω) in 1 0 (Ω)), we can show easily the convergence of the linear terms except the non-linear ones which can be written as: The strong convergence of ℎ to in 1 0 (Ω), and in 6 (Ω), the strong convergence of u ℎ toū in 3 (Ω) , and the weak convergence of ℎ to¯in 1 0 (Ω) lead to the convergence of (3.29).

30)
then, we have the following a priori error estimates:

33)
and (1) Let us estimate the velocity error in terms of the temperature error. By taking the difference between the first equations of ( ) and ( ℎ,1 ) and testing with v = v ℎ ∈ ℎ , we obtain Then, by inserting ℱ ℎ (u) and testing with v ℎ = ℱ ℎ (u) − u ℎ that belongs indeed to ℎ , we easily derive Let us bound the second term in the right hand side of (3.36). We have (3.37) Then, Thus, the monotonicity of and the fact that f 1 is f1 -Lipschitz with values in IR allow us to obtain,

(3.39)
To treat the last inequality, we bound all the terms of the right hand side containing = ‖ v ℎ ‖ 2 (Ω) using the formula and, the term containing =‖ v ℎ ‖ 3 (Ω) , using the formula Then, we infer the following bound: (3.40) By using the following triangle inequality Thus, a triangle inequality allows us to get Hence relations (3.32) and (3.33) are proved (2) The proof of the error estimate for the pressure follows the same lines of the previous step. By taking the difference between the first equations of ( ) and ( ℎ ), inserting ℎ ( ) and testing with v ℎ ∈ ℎ , we obtain (3.44) In order to estimate the second term in the right hand side of (3.44), we will proceed as follows: . By following the same steps of the previous part, and the inf-sup condition (3.3), we get . (3.47) The following triangle inequality allows us to get (3.50) The terms in the last two lines of the right-hand side are bounded by (3.51) Then the choice ℎ = ℎ ( ) − ℎ , the antisymmetric property of the transport term, the fact that u ℎ is bounded in 3 (Ω) as and Sobolev's embedding yield )︂ .
(3.52) By using the following triangle inequality: By using the properties of the operators ℎ , ℱ ℎ and ℎ , we get the following result: Theorem 3.7. Under the assumptions of Theorem 3.6 and if the solution (u, , ) of problem ( ) satisfies ∈ 2 (Ω), u ∈ 1,3 (Ω) and ∈ 2 (Ω), then we have the following a priori error estimates: where 1 and 2 are strictly positive constants independent of ℎ.

Iterative algorithm
In order to solve the discrete system, we propose in this section an iterative algorithm which converges to the exact solution under additional conditions on the exact solution.
The algorithm proceeds as follows: Let u ℎ 0 ∈ ℎ and 0 ℎ ∈ 0 the initial guesses. Having (u ℎ , ℎ ) ∈ ℎ × ℎ at each iteration , we compute (u +1 where is a real strictly positive parameter. Later on, the parameter will be chosen to ensure the convergence of algorithm ( ℎ ). At each iteration , having u ℎ and ℎ , the first two lines of ( ℎ ) computes (u +1 ℎ , +1 ℎ ). Next, we substitute u +1 ℎ by its value in the third equation of ( ℎ ) to compute +1 ℎ .
In the following, we study Scheme ( ℎ ), and we begin by proving the existence and uniqueness of the corresponding solution.
The bound (4.2) can be deduced immediately by taking ℎ = +1 ℎ in the third equation of problem ( ℎ ), and by using the Cauchy-Schwartz inequality.
To prove the bound (4.6), we need first to estimate the error By inserting u ℎ in the second and third terms of the last equation, we get, (4.9) Using the properties of −1 , the Cauchy-Schwartz inequality and relation (3.2) give the following (4.10) We apply the relation ≤ 1 2 2 + 2 2 with = 3 to each term on the right-hand side of the previous inequality, and we obtain (4.11) Therefore, Assumption 2.1 and Relation (4.2) allow us to get where (4.13) Then, we are now in position to show relation (4.6). We consider the first equation of problem ( ℎ ) with v ℎ = u +1 ℎ , and we obtain ∫︁ (4.14) Using the properties of −1 , the Cauchy-Shwartz inequality and the relations ≤ 1 2 2 + 2 2 and 2 ≤ )︁ wih = and = (︁ 3 4

(4.15)
We denote by )︀ and we get by using (4.12) that Therefore, we obtain the following bound: (4.16) We now prove estimate (4.6) by induction on ≥ 1 under some condition on . Starting with relation (4.3), we suppose that we have We have two situations: . By using the induction condition (4.17), taking 2 > 16 )︀ , )︀ > 0, and deduce from relation (4.16) that Then relation (4.6) holds. The bound (4.7) is a simple consequence of (4.16) and (4.6).
The next theorem shows the convergence of the solution (u ℎ , ℎ , ℎ ) of problem ( ℎ ) to the solution of problem ( ℎ ).
Theorem 4.2. In addition to Assumption 2.1, we assume that the concentration solution of the problem ( ) satisfies Under the assumptions of Theorem 4.1, and if satisfies the condition where 1 is the constant in (3.54), then the solution (u ℎ , ℎ , ℎ ) of problem ( ℎ ) converges in 2 (Ω) × 2 (Ω) × 1 (Ω) to the solution of problem ( ℎ ).
Proof. We start by subtracting the third equation of problem ( ℎ ) from the one of problem ( ℎ ) to get Inserting ∇ ℎ in the last term of the left-hand side and ℎ in the first term of the right-hand side of the previous relation lead to (4.23) Finally, by inserting ∇ in the second term of right-hand side and using the Green's formula and the antisymmetric property, we get (4.24) By taking ℎ = ℎ − +1 ℎ , we obtain and then we deduce the convergence of the sequence (u +1 ℎ − u ℎ ) in 2 (Ω) , and then the convergence of the sequence u ℎ in 2 (Ω) . By taking the limit of (4.32), we get that u +1 ℎ converges to u ℎ in 2 (Ω) . Relation (4.26) allows us to deduce that +1 ℎ converges to ℎ in 1 0 (Ω). Next, we study the convergence of the pressure, taking the difference between the first equations of systems ( ℎ ) and ( ℎ ), we obtain for all v ℎ ∈ ℎ the equation We get by using the inverse inequality (3.2) the following: For a given mesh, owning the inf-sup condition (3.3), and using the strong convergence of u ℎ to u ℎ in 2 (Ω) , and of ℎ to ℎ in 1 0 (Ω), we deduce the strong convergence of ∇ ℎ to ∇ ℎ in 3 2 (Ω). Furthermore, the fact that ℎ and ℎ are in the discrete space of IP 1 finite elements ℎ ⊂ 2 0 (Ω) which is defined in (3.6), allows us to deduce the strong convergence of ℎ to ℎ in 2 (Ω).

Remark 4.3.
If ℎ is not small enough, we can replace the conditions (4.19) and (4.21) in Theorem 4.2 by the following condition: In fact, this condition can be used in the relation (4.31) in the proof of the previous theorem.  .4) and (4.20) satisfied by to ensure the convergence of the iterative solution (u ℎ , ℎ , ℎ ) to the discrete solution (u ℎ , ℎ , ℎ ) lead us to consider as a function of ℎ ( (ℎ)) and satisfying these two relations. Thanks to the triangle inequality, we have: Under the assumptions of Theorem 4.2 with (ℎ) satisfying conditions (4.4) and (4.20), (u ℎ , ℎ , ℎ ) converges to (u ℎ , ℎ , ℎ ) in ℎ × ℎ × ℎ , and, there exists an integer 0 (ℎ) depending on ℎ such that for all ≥ 0 (ℎ) we have Consequently, for all ≥ 0 (ℎ), we get Thus, we obtain the convergence of the iterative solution to the exact one.

Numerical results
In this section, we present numerical experiments corresponding to our coupled problem for = 2. These simulations have been performed using the code FreeFem++ due to Hecht and Pironneau (see [25]).
We will show in this section numerical investigations corresponding to problems ( ℎ ) by using for the convergence the stopping criterion Err ≤ where is a given tolerance considered in this work equal to 10 −5 and Err is defined by The initial guesses u 0 ℎ and 0 ℎ are considered in one of these two situations: (1) 0 ℎ = 0 and u 0 ℎ = 0. (2) 0 ℎ = 0 and u 0 ℎ = u 0 ℎ are calculated by using Darcy's problem which corresponds to = = 0. We will see later that the second case where u 0 ℎ is the solution of Darcy's problem improve the convergence of the algorithms.
We also consider the errors which describe the rate of the a priori error estimation for a large values of the iteration index .

First numerical test: analytical solution
In this section, we will show numerical results corresponding to the problem where we know the exact solution. Let Ω =]0, 1[ 2 ⊂ IR 2 where each edge is divided into equal segments so that Ω is divided into 2 equal squares and finally into 2 2 equal triangles. For simplicity, we take = = 1.

Case where =
In this part, we take = .
To study the dependency of the convergence with the parameter , we consider = 60, = 20, = 20, = 0 = 1, and for each , we stop the algorithm ( ℎ ) when the error Err < 1 −5 . Tables 1 shows the error (Err) and the number of iterations (Nbr) for u 0 ℎ = 0 and 0 ℎ = 0, while Table 2 shows similar results for u 0 = u 0 ℎ and 0 ℎ = 0. These two tables describe the convergence of algorithm ( ℎ ) with respect to . We remark that the number of iterations is relatively small when is large. In both cases, the best convergence is obtained for = 100. The main advantage is for the case where u 0 ℎ = u 0 ℎ are computed with Darcy's problem is that the number of iterations (Nbr) is less than the one obtained with the case u 0 ℎ = 0. For further studies, we consider = 10, = 20, = 100 and the initial guesses u 0 ℎ = u 0 ℎ and 0 ℎ = 0. Tables 3 and 4 show the obtained rate of convergence which seems in agreement with the theoretical findings. We notice that the theoretical rate of convergence of the velocity in norm 2 (Ω) 2 , the pressure in norm 1, 3 2 (Ω) and the concentration in norm 1 0 (Ω) are equal to 1; the rate of convergence of the velocity in norm 3 (Ω) 3 is 2/3.  Table 3. Rate of convergence of the velocity in norms 2 (Ω) 2 and 3 (Ω) 3 . Example (5.1) with algorithm ( ℎ ) ( = 20, = 20).   3 + sin( ) sin( ) )︂ .
We follow the same numerical tests performed above and we consider = 60, = 20, = 20, = 0 = 1. Based on the previous test case, we will perform here numerical simulations with u 0 ℎ = u 0 ℎ and 0 ℎ = 0. Table  5 describes the convergence of algorithm ( ℎ ) with respect to . In this case, the best convergence is also obtained for = 100. Let us now study the rate of convergence of the errors. we consider = 10, = 20, = 100 and the initial guesses u 0 ℎ = u 0 ℎ and 0 ℎ = 0. Tables 6 and 7 show that the numerical rate of convergence seems in agreement with the theoretical findings.

Comparison between Darcy and Darcy-Forchheimer
In this section, we treat numerical test (taken from [32]) showing the difference between Darcy and Darcy-Forchheimer systems. We take = 60, = 10, In the following, we will compare the numerical velocity, pressure and concentration corresponding to Darcy-Forchheimer Problem (for = 10) and Darcy Problem (for = 0). Figures 1-6 show that there are differences between the distributions and values of the numerical velocities, pressures and concentration. The biggest difference between the Darcy and the Darcy-Forchheimer is in the values of the velocity (which was expected).

Second numerical test: Driven cavity
The driven cavity is a standard benchmark for testing the performance of algorithms in fluid problems. It is treated in several works (see [8,16,27,34]). In this section, we show numerical simulation corresponding to the Lid Driven Cavity in order to study the dependency of the convergence with respect to and the data.
Let Ω =]0, 1[ 2 , = , = 1, 0 = 0, = 20, f 0 = 0, f 1 ( ) = (10 , 10 ), and = 0. We complete the Darcy-Forchheimer equations with the boundary conditions u.n = 0 on Ω, and the concentration equation with the boundary condition = ( is a parameter) on Γ 1 = [0, 1] × {1} (top of Ω), and = 0 on Ω∖Γ 1 . In this section, the initial guesses of algorithm ( ℎ ) are 0 ℎ = 1 and u 0 ℎ = u 0 ℎ . We begin first by testing the convergence of the algorithm with respect of for a given = 20. We consider = 20 and we test the algorithm for multiple values of . We consider that the algorithm does not converge  if the condition Err < 1 −5 is not reached after 5000 iterations. Table 8 shows for = 20, the dependency of the convergence of the algorithm with respect to and the better convergence corresponds to the value = 10. This result shows clearly that the convergence depends on as announced in relation (4.20) of Theorem 4.2. Figures 7-9 show the velocity, pressure, and concentration in Ω for = 10 and = 20.
Let us now test the convergence with respect to for = 10. Table 9 shows the dependency of the convergence of the algorithm with respect to (i.e. with respect to the concentration ). This result shows clearly that the convergence depends on the concentration as announced in relation (4.19) of Theorem 4.2.          1  20  100  150  170  175  180 200   Nbr  24  27  77  235 728 1493 -conv/div conv conv conv conv cov conv div div