NEW MIXED FINITE ELEMENT METHODS FOR THE COUPLED STOKES AND POISSON–NERNST–PLANCK EQUATIONS IN BANACH SPACES

. In this paper we employ a Banach spaces-based framework to introduce and analyze new mixed finite element methods for the numerical solution of the coupled Stokes and Poisson–Nernst– Planck equations, which is a nonlinear model describing the dynamics of electrically charged incompressible fluids. The pressure of the fluid is eliminated from the system (though computed afterwards via a postprocessing formula) thanks to the incompressibility condition and the incorporation of the fluid pseudostress as an auxiliary unknown. In turn, besides the electrostatic potential and the concentration of ionized particles, we use the electric field (rescaled gradient of the potential) and total ionic fluxes as new unknowns. The resulting fully mixed variational formulation in Banach spaces can be written as a coupled system consisting of two saddle-point problems, each one with nonlinear source terms depending on the remaining unknowns, and a perturbed saddle-point problem with linear source terms, which is in turn additionally perturbed by a bilinear form. The well-posedness of the continuous formulation is a consequence of a fixed-point strategy in combination with the Banach theorem, the Babuˇska–Brezzi theory, the solvability of abstract perturbed saddle-point problems, and the Banach– Neˇcas–Babuˇska theorem. For this we also employ smallness assumptions on the data. An analogous approach, but using now both the Brouwer and Banach theorems, and invoking suitable stability conditions on arbitrary finite element subspaces, is employed to conclude the existence and uniqueness of solution for the associated Galerkin scheme. A priori error estimates are derived, and examples of discrete spaces that fit the theory, include, e.g. , Raviart–Thomas elements of order 𝑘 along with piecewise polynomials of degree ď 𝑘 . In addition, the latter yield approximate local conservation of momentum for all three equations involved. Finally, rates of convergence are specified and several numerical experiments confirm the theoretical error bounds. These tests also illustrate the aforementioned balance-preserving properties and the applicability of the proposed family of methods.


Introduction
Fluid mixtures with electrically charged ions are critical for many industrial processes and natural phenomena.Notable examples of current interest are efficient energy storage and electrodialysis cells, design of nanopore sensors, electro-osmotic water purification techniques, and even drug delivery in biological tissues [40].One of the most well-known models for liquid electrolytes is the Poisson-Nernst-Planck/Stokes system.It describes the isothermal dynamics of the molar concentration of a number of charged species within a solvent.This classical model is valid for the regime of relatively small Reynolds numbers and it is written in terms of the concentrations, the barycentric velocity of the mixture, the pressure of the mixture, and the electrostatic potential.The system is strongly coupled and the set of equations consist of the transport equations for each dilute component of the electrolyte, a diffusion equation for the electrostatic equilibrium, the momentum balance for the mixture (including a force exerted by the electric field), and mass conservation.
Solving these systems lends itself difficult due to coupling nonlinearities of different nature.Numerical methods for incompressible flow equations coupled with Poisson-Nernst-Planck equations that are based on finite element schemes in primal formulation (also including stabilized and goal-adaptive methods) can be found in [3,20,32,35,37,39], finite differences in e.g., [34], finite volume schemes in [38], spectral elements in [36], and also for virtual element methods in [17].Regarding formulations using mixed methods, the first works addressing Stokes/PNP systems are relatively recent [27,28], where a stabilized mixed method is employed for the Poisson problem, whereas the usual primal approach is applied to the Stokes, Navier-Stokes and Nernst-Planck equations, and all them within a Hilbertian framework.Mixed variational formulations are particularly interesting when direct discrete approximations of further variables of physical relevance are required.A recent approach to mixed methods consists in defining the corresponding variational settings in terms of Banach spaces instead of the usual Hilbertian framework, and without augmentation.As a consequence, the unknowns belong now to the natural spaces that are originated after carrying out the respective testing and integration by parts procedures, simpler and closer to the original physical model formulations arise, momentum conservative schemes can be obtained, and even other unknowns can be computed by postprocessing formulae.As a non-exhaustive list of contributions taking advantage of the use of Banach frameworks for solving the aforementioned kind of problems, we refer to [4,[8][9][10][11]13,14,24,26,29], and among the different models considered there, we find Poisson, Brinkman-Forchheimer, Darcy-Forchheimer, Navier-Stokes, chemotaxis/Navier-Stokes, Boussinesq, coupled flow-transport, and fluidized beds.Nevertheless, and up to our knowledge, no mixed methods with the described advantages seem to have been developed so far for the coupled Stokes and Poisson-Nernst-Planck equations.
As motivated by the previous discussion, the goal of this paper is to develop a Banach spaces-based formulation yielding new mixed finite element methods for, precisely, the coupled Stokes and Poisson-Nernst-Planck equations.The main novelties with respect to [27,28] refer to the use of mixed methods for each one of the equations involved, the setting of the resulting variational formulation within a Banach framework, and the no need of incorporating any additional stabilization term.The rest of the manuscript is organized as follows.Required notations and basic definitions are collected at the end of this introductory section.In Section 2 we describe the model of interest and introduce the additional variables to be employed.The mixed variational formulation is deduced in Section 3.After some preliminaries, the respective analysis is split according to the three equations forming the whole system.In particular, the right spaces to which the trial and test functions must belong are derived in each case by applying suitable integration by parts formulae jointly with the Cauchy-Schwarz and Hölder inequalities.In Section 4 we utilize a fixed-point approach to study the solvability of the continuous formulation.The Babuška-Brezzi theory and recent results on perturbed saddle-point problems, both in Banach spaces, along with the Banach-Nečas-Babuška theorem, are utilized to prove that the corresponding uncoupled problems are well-posed.The classical Banach fixed-point theorem is then applied to conclude the existence of a unique solution.In Section 5 we proceed analogously to Section 4 and, under suitable stability assumptions on the discrete spaces employed, show existence and then uniqueness of solution for the Galerkin scheme by applying the Brouwer and Banach theorems, respectively.A priori error estimates are also derived here.Next, in Section 6 we define explicit finite element subspaces satisfying those conditions, and provide the associated rates of convergence.Finally, several numerical examples confirming the latter, showing the good performance of the method, and illustrating the approximate local conservation of momentum, are reported in Section 7.

Preliminary notations
Throughout the paper, Ω is a bounded Lipschitz-continuous domain of   ,  P t2, 3u, which is star shaped with respect to a ball, and whose outward normal at Γ :" BΩ is denoted by .Standard notation will be adopted for Lebesgue spaces   pΩq and Sobolev spaces  , pΩq and  , 0 pΩq, with  ě 0 and  P r1, `8q, whose corresponding norms, either for the scalar and vectorial case, are denoted by }¨} 0,;Ω and }¨} ,;Ω , respectively.Note that  0, pΩq "   pΩq, and if  " 2 we write   pΩq instead of  ,2 pΩq, with the corresponding norm and seminorm denoted by }¨} ,Ω and | ¨|,Ω , respectively.In addition, letting ,  1 P p1, `8q conjugate to each other, that is such that 1{ `1{ 1 " 1, we denote by  1{ 1 , pΓq the trace space of  1, pΩq, and let  ´1{ 1 , 1  pΓq be the dual of  1{ 1 , pΓq endowed with the norms }¨} ´1{ 1 , 1 ;Γ and }¨} 1{ 1 ,;Γ , respectively.On the other hand, given any generic scalar functional space  , we let M and M be the corresponding vectorial and tensorial counterparts, whereas }¨} will be employed for the norm of any element or operator whenever there is no confusion about the spaces to which they belong.Furthermore, as usual, I stands for the identity tensor in R :"  ˆ , and | ¨| denotes the Euclidean norm in   .Also, for any vector field v " p  q "1, we set the gradient and divergence operators, respectively, as ∇v :" ´B B ¯,"1, and divpvq :" Additionally, for any tensor fields  " p  q ,"1, and  " p  q ,"1, , we let divp q be the divergence operator div acting along the rows of  , and define the transpose, the trace, the tensor inner product operators, and the deviatoric tensor, respectively, as  t " p  q ,"1, , trp q "  ÿ "1   ,  :  :"  ÿ ,"1     , and  d :"  ´1  trp qI.

The model problem
We consider the nonlinear system given by the coupled Stokes and Poisson-Nernst-Planck equations, which constitute an electrohydrodinamic model describing the stationary flow of a Newtonian and incompressible fluid occupying the domain Ω Ď   ,  P t2, 3u, with polygonal (resp.polyhedral) boundary Γ in  2 (resp. 3 ).Under the assumption of isothermal properties, equal molar volumes and molar masses for each species, the behavior of the system is determined by the concentrations  1 and  2 of ionized particles, and by the electric current field .Mathematically speaking, and firstly regarding the fluid, we look for the barycentric velocity u and the pressure  of the mixture, such that pu, q is solution to the Stokes equations ´∆u `∇ " ´p 1 ´2 q ´1 `f in Ω, divpuq " 0 in Ω, u " g on Γ, where  ą 0 is the constant viscosity,  is the heterogeneous dielectric coefficient, also known as the electric conductivity coefficient, f is a source term, g is the Dirichlet datum for u on Γ, and the null mean value of  has been incorporated as a uniqueness condition for this unknown.In addition, ,  1 and  2 solve the Poisson-Nernst-Planck equations, which depend on the velocity u and are given by  " ∇ in Ω, ´divpq " p 1 ´2 q ` in Ω, where  is the electrostatic potential, and for each  P t1, 2u where  1 and  2 are the diffusion coefficients,   :" We stress that in order to solve (3), u and  are needed.In turn, (1) requires  1 ,  2 and , whereas (2) makes use of  1 and  2 .This multiple coupling is illustrated through the graph provided in Figure 1, where the vertexes represent the aforementioned equations and the arrows, properly labeled with the unknowns involved, show the respective dependence relationships.
Furthermore, since we are interested in employing a fully mixed variational formulation for the coupled model ( 1)-(3), we introduce the auxiliary variables of pseudostress and, for each  P t1, 2u, the total (diffusive, cross-diffusive, and advective) ionic fluxes Thus, applying the matrix trace in (5) and using the incompressibility condition, we deduce that  " ´1  trpq, (7) so that, incorporating the latter expression into (5),  is eliminated from the system (1)-( 3), which can then be rewritten in terms of the unknowns , u, , ,   and   ,  P t1, 2u, as We notice here that the uniqueness condition for  has been rewritten equivalently as the null mean value constraint for trpq.

The fully mixed formulation
In this section we derive a suitable Banach spaces-based variational formulation for (8) by splitting the analysis in four sections.The first one collects some preliminary discussions and known results, and the remaining three deal with each one of the pairs of equations forming the whole nonlinear coupled system (8), namely Stokes, Poisson, and Nernst-Planck.

Preliminaries
We begin by noticing that there are three key expressions in (8) that need to be looked at carefully before determining the right spaces where the unknowns must be sought, namely p 1 ´2 q  ´1,      ´1 and  ´1    u in the first and fifth rows of (8).More precisely, ignoring the bounded above and below functions  ´1 and  ´1  , as well as the constant   , and given test functions v and   associated with u and   , respectively, straightforward applications of the Cauchy-Schwarz and Hölder inequalities yield and similarly ˇˇˇż where ℓ,  P p1, `8q are conjugate to each other.In this way, denoting ´1 (conjugate of ),  :" 2, and  :" it follows that the above expressions make sense for   P   pΩq, , u P L  pΩq, and v,   P L 2 pΩq.The specific choice of ℓ, and hence of , ,  and the respective conjugates  and , will be addressed later on, so that meanwhile we consider generic values for the indexes defined in (10).Having set the above preliminary choice for the space to which  belongs, we deduce from the first equation in the third row of (8) that  should be initially sought in  1, pΩq.In turn, using that  1 pΩq is embedded in   pΩq for  P r1, `8q in  2 (resp. P r1, 6s in  3 ), and for reasons that will become clear below, the unknowns   ,  P t1, 2u, and u are initially sought in  1 pΩq and H 1 pΩq, respectively, certainly assuming that  and  verify the indicated ranges, namely ,  P p2, `8q in  2 , and ,  P p2, 6s in  3 .Note that in terms of ℓ the latter constraint becomes ℓ P r 3  2 , 3s, which yields  P r3, 6s.Equivalently,  P r 3 2 , 3s and  P r3, 6s, though going through the respective intervals in the opposite direction to ℓ and , respectively.

The Stokes equations
Let us first notice that, applying (14) with  "  to  P Hpdiv  ; Ωq and u P H 1 pΩq, and using the Dirichlet boundary condition on u, for which we assume from now on that g P H 1{2 pΩq, we obtain ż Ω  : ∇u " ´żΩ u ¨divp q `x , gy, and thus, the testing of the first equation in the first row of (8) against  yields Note from the second term on the left-hand side of ( 16) that, knowing that divp q P L  pΩq, it actually suffices to look for u in L  pΩq, which is coherent with a previous discussion on the space to which this unknown should belong.In addition, testing the second equation in the first row of (8) against v P L  pΩq, for which we require that f P L  pΩq, we get ż which makes sense for divpq P L  pΩq.Hence, due to the last equation in the second row of (8), it follows that we should look for  in H 0 pdiv  ; Ωq, where H 0 pdiv  ; Ωq :" "  P Hpdiv  ; Ωq : Moreover, it is easily seen that there holds the decomposition and that the incompressibility of the fluid forces the compatibility condition on g given by As a consequence of the above, we realize that imposing (16) for each  P Hpdiv  ; Ωq is equivalent to doing it for each  P H 0 pdiv  ; Ωq.Furthermore, since  ą 2 it follows that L  pΩq is embedded in L 2 pΩq, which, along with the estimate (9a), confirms that the first term on the right-hand side of ( 17) is also well-defined.In this way, denoting from now on  :" p 1 ,  2 q, and joining ( 16) and (17), we arrive at the following mixed variational formulation for the Stokes equations (given by the first two rows of ( 8)): Find p, uq P H ˆQ such that ap,  q `bp , uq " Fp q @ P H, bp, vq " G , pvq @v P Q, where H :" H 0 pdiv  ; Ωq, and the bilinear forms a : H ˆH Ñ  and b : H ˆQ Ñ , and the functional F : H ÝÑ , are defined, respectively, as bp , vq :" Fp q :" x , gy whereas, given  :" p 1 ,  2 q P L  pΩq and  P L  pΩq, the functional G , : Q ÝÑ  is set as It is readily seen that, endowing H with the corresponding norm from (12b), that is and recalling that } ¨}0,;Ω is that of Q, the bilinear forms a and b, and the linear functionals F and  , , are all bounded.Indeed, applying the Cauchy-Schwarz and Hölder inequalities, noting that } d } 0,Ω ď } } 0,Ω for all  P H, invoking the identity ( 14) along with the continuous injection i  : H 1 pΩq Ñ L  pΩq, using (9a) together with the fact that } ¨}0,Ω ď |Ω| p´2q{2 } ¨}0,;Ω , and bounding  ´1 according to (4), we deduce the existence of positive constants, denoted and given as }a} :" 1  , }b} :" 1, }F} :" p1 `}i  }q }g} 1{2,Γ , and }G} :" max We end this section by emphasizing, according to the previous discussion, that the introduction of the pseudostress  as an auxiliary unknown leads to the derivation of simple postprocessing formulae for the pressure  (cf.(7)) and the velocity gradient ∇u (cf.first equation in the first row of ( 8)).In addition, it allows us to seek the velocity u in a Lebesgue space, which is certainly less regular, whence the corresponding finite element subspace, not requiring any continuity property, can be chosen cheaper and easier to implement.

The electrostatic potential equations
We begin the derivation of the mixed formulation for the Poisson equation by testing the first equation in the third row of (8) against  P H  pdiv  ; Ωq.In this way, applying (15) with  "  and  1 "  to the given  and  P  1, pΩq, and employing the Dirichlet boundary condition on , for which we assume that  P  1{, pΓq, we get ż In turn, testing the second equation in the third row of (8) against  P   pΩq, which requires to assume that  P   pΩq, we obtain ż which certainly makes sense for divpq P   pΩq.Thus, recalling from (9a) and (9b) that  must belong to L  pΩq, it follows from the above that this unknown should be sought then in H  pdiv  ; Ωq.Furthermore, bearing in mind from (9a) to (9c) that  1 and  2 must belong to   pΩq, we notice that in order for the first term on the right-hand side of (26) to make sense, we require that  ě , which is assumed from now on.Therefore, placing together (25) and (26), we obtain the following mixed variational formulation for the electrostatic potential equations (given by the third and fourth rows of ( 8)): Find p, q P  2 ˆ1 such that p, q `1 p, q "  pq @ P  1 ,  2 p, q "   pq @ P  2 , where  2 :" H  pdiv  ; Ωq,  1 :"   pΩq,  1 :" H  pdiv  ; Ωq,  2 :"   pΩq, and the bilinear forms  :  2 ˆ1 Ñ  and   :   ˆ Ñ ,  P t1, 2u, and the functional  :  1 Ñ , are defined, respectively, as p, q :" p, q :" pq :" x ¨, y Γ @ P  1 , (29c) whereas, given  :" p 1 ,  2 q P L  pΩq, the functional   :  2 Ñ  is defined by We end this section by establishing the boundedness of ,   ,  P t1, 2u,  , and   , for which we recall that the norms of  1 and  2 are defined by (12c) with  "  and  " , respectively, whereas those of  1 and  2 are certainly given by } ¨}0,;Ω and } ¨}0,;Ω , respectively.Then, employing again the Cauchy-Schwarz and Hölder inequalities, bounding  ´1 according to (4), and using that } ¨}0,;Ω ď |Ω| p´q{ } ¨}0,;Ω , which follows from the fact that  ě , we find that there exist positive constants }} :" 1  0 , } 1 } " } 2 } :" 1, and }} :" max such that |p, q| ď }} }} 2 }} 1 @p, q P  2 ˆ1 , |  p, q| ď }  } }}  }}  @p, q P   ˆ1 , @ P t1, 2u, and Regarding the boundedness of  , we need to apply Lemma A.36 of [19], which, along with the surjectivity of the trace operator mapping  1, pΩq onto  1{, pΓq, yields the existence of a fixed positive constant   , such that for the given  P  1{, pΓq, there exists   P  Hence, employing (15) with p,  1 q " p, q and p , q " p,   q, and then using Hölder's inequality, we arrive at with } } :"   }} 1{,;Γ .

The ionized particles concentration equations
We now deal with the Nernst-Planck equations, that is the fifth and sixth rows of ( 8), for which we proceed analogously as we did for the Stokes equations in Section 3.2.More precisely, applying (13) with  "  to   P Hpdiv  ; Ωq and   P  1 pΩq, and using the Dirichlet boundary condition on   , for which we assume from now on that   P  1{2 pΓq, we obtain so that the testing of the equation in the fifth row of (8) against   , yields Since divp  q P   pΩq, we notice from the second term on the left-hand side of (34) that it suffices to look for   in   pΩq, which, similarly as for Stokes, is coherent with a previous discussion on where to seek this unknown.
Notice here that the second, fourth, and sixth rows of (41) constitute the conservation of momentum for each respective equation.We will refer again to this subject from the discrete point of view later on in Sections 5.1 and 6.2.
We end this section by stressing that, as compared with previously studied Banach spaces-based mixed formulations for other coupled nonlinear models (see, e.g., [4,14,24,26]), the main novelty of the analysis to be developed for (41) has to do with the occurrence in its last two rows of the perturbed saddle point scheme in Banach spaces represented by the bilinear forms   ,   , and   .Indeed, up to our knowledge, the present one constitutes the first work applying the theoretical results provided recently in [15] to perform the continuous and discrete analyses of a problem showing that structure.

The continuous solvability analysis
In this section we proceed as in several related previous contributions (see, e.g., [11] and the references therein), and employ a fixed-point strategy to address the solvability of (41).

The fixed-point strategy
In order to rewrite (41) as an equivalent fixed point equation, we introduce suitable operators associated with each one of the three problems forming the whole nonlinear coupled system.Indeed, we first let p  : p 1 ˆ2 q ˆ2 Ñ Q be the operator defined by p  p, q :" p u @p, q P p 1 ˆ2 q ˆ2 , where pp , p uq P H ˆQ is the unique solution (to be confirmed below) of problem (19) (equivalently, the first two rows of (41)) with p, q instead of p, q, that is app ,  q `bp , p uq " Fp q @ P H, bpp , vq " G , pvq @v P Q.
In turn, we let s  :  1 ˆ2 Ñ  2 be the operator given by s  pq :" s  @ P  1 ˆ2 , where p s , s q P  2 ˆ1 is the unique solution (to be confirmed below) of problem ( 27) (equivalently, the third and fourth rows of (41)) with  instead of , that is p s , q `1 p, s q "  pq @ P  1 ,  2 p s , q "   pq @ P  2 . (43) Similarly, for each  P t1, 2u, we let r   :  2 ˆQ Ñ   be the operator defined by where pr   , r   q P   ˆ is the unique solution (to be confirmed below) of problem (36) (equivalently, the fifth and sixth rows of (41)) with p, vq instead of p, uq, that is so that we can define the operator r  :  2 ˆQ Ñ p 1 ˆ2 q as: Finally, defining the operator T : we observe that solving (41) is equivalent to seeking a fixed point of T, that is: Find  P  1 ˆ2 such that Tpq " . (47)

Well-definedness of the operator p 𝑇
Here we apply Theorem 2.34 of [19] to show that, given an arbitrary p, q P `1 ˆ2 ˘ˆ 2 , (42) is wellposed, equivalently that p  is well-defined.We remark that p, q only influences the functional G , , and that the boundedness of all the bilinear forms and linear functionals defining (42), has already been established in (23) and (24).Hence, the discussion below just refers to the remaining hypotheses to be satisfied by a and b.We begin by letting V be the kernel of the operator induced by b, that is V :" t P H : bp , vq " 0 @v P Qu, which, according to the definitions of H, Q, and b (cf.( 20), (21b)), along with the fact that L  pΩq is isomorphic to the dual of L  pΩq, yields V :" t P H 0 pdiv  ; Ωq : divp q " 0u.
Next, we recall that a slight modification of the proof of Lemma 2.3 from [22] allows to prove that for each  ě 2 `2 (see, e.g., [8], Lem.3.1 for the case  " 4{3, which is extensible almost verbatim for any  in the indicated range) there exists a constant   , depending only on Ω, such that Then, assuming that  ě 2 `2 , and using (49), we deduce from the definition of a (cf.(21a)), and similarly to Lemma 3.2 of [8], that ap ,  q ě  } } 2 div;Ω with  :"   {.Hence, thanks to (50), it is straightforward to see that a satisfies the hypotheses specified in Theorem 2.34, equation (2.28) of [19] with the foregoing constant .In order to fulfill all the hypotheses of the latter theorem, and knowing from ( 23) and ( 24) that the boundedness of the corresponding bilinear forms and linear functionals has already been established, it only remains to show the continuous inf-sup condition for b.Moreover, being this result already proved for the particular case  " 4{3 (cf.[8], Lem.3.3 and [24], Lem.3.5 for a closely related one), and arising no significant differences for an arbitrary  ě 2 `2 , we provide below, and for sake of completeness, only the main aspects of its proof.
Indeed, given v P Q :" L  pΩq, we first recall from (10) that  ą 2, and set v  :" |v| ´2 v, which is easily seen to satisfy In what follows, we make use of both, the Poincaré inequality, which refers to the existence of a positive constant  P , depending on Ω, such that  P }z} 2  1,Ω ď |z| 2 1,Ω @z P H 1 0 pΩq, and the continuous injection i  : H 1 pΩq Ñ L  pΩq for the indicated range of .Then, we let w P H 1 0 pΩq be the unique solution of: ş Ω ∇w ¨∇z " ´şΩ v  ¨z for all z P H 1 0 pΩq, which is guaranteed by the classical Lax-Milgram Lemma, and notice, thanks to the corresponding continuous dependence estimate, that }w} 1,Ω ď }i} P }v  } 0,;Ω .Hence, defining  :" ∇w P L 2 pΩq, we deduce that divpq " v  in Ω, so that  P Hpdiv  ; Ωq, and }} div;Ω ď `1 `}i} P ˘}v  } 0,;Ω .Finally, letting  0 be the H 0 pdiv  ; Ωq-component of , it is clear that divp 0 q " v  and that } 0 } div;Ω ď }} div;Ω , whence bounding by below with  :"  0 P H, and using the definition of b (cf.(21b)) along with the above identities and estimates, we conclude that sup with  :" `1`} i} P ˘´1 .The foregoing inequality (51) proves Theorem 2.34, equation (2.29) of [19] and completes the hypotheses of this theorem.
Consequently, the well-definedness of the operator p  is stated as follows.
Theorem 4.1.For each p, q P p 1 ˆ2 q ˆ2 there exists a unique pp , p uq P H ˆQ solution to (42), and hence one can define p  p, q :" p u P Q.Moreover, there exists a positive constant  p  , depending only on , }i  },  0 , |Ω|, , and , and hence independent of p, q, such that Proof.Given p, q P p 1 ˆ2 q ˆ2 , a direct application of Theorem 2.34 from [19] guarantees the existence of a unique pp , p uq P H ˆQ solution to (42).Then, the corresponding a priori estimate in Theorem 2.34, equation (2.30) of [19] gives which, according to the identities and estimates given by ( 23) and ( 24), along with some algebraic manipulations, yields (52) and finishes the proof.
Regarding the a priori bound for the component p  of the unique solution to (42), it follows from Theorem 2.34, equation (2.30) of [19] that which yields the same inequality as ( 52), but with a different constant.Hence, choosing the largest of the respective constants, and still denoting it by  p  , we can summarize the a priori estimates for p u and p  by saying that both are given by the right-hand side of (52).

Well-definedness of the operator s 𝑇
We now employ Theorem 2.1, Section 2.1 of [5] to prove that, given an arbitrary  P  1 ˆ2 , ( 43) is well-posed, equivalently that s  is well-defined.Similarly as for Section 4.2.1, we first stress that  is utilized only to define the functional   , and that the boundedness of all the bilinear forms and functionals defining (43), was already established by (30) and (31).In this way, it only remains to show that ,  1 , and  2 satisfy the corresponding hypotheses from Theorem 2.1, Section 2.1 of [5].To this end, and because of the evident similarities, we follow very closely the analysis in Section 3.2.3 from [11], which, in turn, suitably adopts the approach from Section 2.4.2 of [26].Indeed, we begin by letting   be the kernel of the operator induced by the bilinear form   , for each  P t1, 2u, that is which, according to the definitions of   and   (cf.( 28)), and   (cf.(29b)), along again with the fact that L  pΩq and L  pΩq can be isomorphically identified with `L pΩq ˘1 and `L pΩq ˘1, respectively, gives and  2 :" t P H  pdiv  ; Ωq : divpq " 0 in Ωu.
Next, in order to establish the inf-sup conditions required for the bilinear form  (cf.[5], Eqs.(2.8) and (2.9)), we resort to Lemma 3.3 of [11], which is recalled below.Lemma 4.2.Let Ω be a bounded Lipschitz-continuous domain of   ,  P t2, 3u, and let ,  1 P p1, `8q conjugate to each other with  (and hence  1 ) lying in " r4{3, 4s if  " 2 r3{2, 3s if  " 3 .Then, there exists a linear and bounded operator In addition, for each z P L  1 pΩq such that divpzq " 0 in Ω, there holds Proof.It reduces to a minor modification of the proof of Lemma 2.3 from [26], for which one needs to apply the well-posedness in  1, pΩq of a Poisson problem with homogeneous Dirichlet boundary conditions (see [23], Thm.3.2 or [30], Thms.1.1 and 1.3 for the vector version of it).The specified ranges for  and  1 are precisely forced by the latter result.We omit further details and refer to the proof of Lemma 3.3 from [11].
Before continuing with the continuous inf-sup conditions for the bilinear forms   ,  P t1, 2u, we now check the feasibility of the indexes employed so far, according to the different constraints that have arisen along the analysis.In fact, from the preliminary discussion provided in Section 3.
In particular, the only possibility for the 3D case is obtained by intersecting the range for  specified in the second row of (63), that is  P r3, 4s, with the one required by Lemma 4.2, that is  P r3{2, 3s, which certainly yields  " 3. The respective conjugate becomes  " 3{2, which clearly verifies  ě 2 `2 " 6{5.The occurrence of this unique way of choosing the exponents does not seem in principle to be coincidental since it has also arose in some related papers when a technical result like Lemma 4.2 (or a similar one), is employed (see, e.g., [26], Eq. (2.20), [11], Eqs.(2.25) and (2.26), and [16], Sect.4.2).However, this is not the case for the stressassisted diffusion problem studied in [25], where the feasible ranges obtained in 3D are actually intervals (see [25], Eqs.(3.70) and (3.71)), and hence it is not possible to conclude a corresponding general rule.
Note that in (64) we have included the consequent ranges for  :"  ´1 and  :"  ´1 as well.However, we remark that the above indexes are not chosen independently, but once  (or its conjugate ) is chosen, then all the remaining ones are fixed.In this regard, and extending a related comment made in Section 3.1, we stress here that in the 2D case the values of the feasible exponents  and  (equivalently, the indexes  and ℓ) vary in opposite directions, namely as  increases,  decreases, and conversely.Similarly, being p, q and p, q conjugate pairs, as the first component of each increases, the second one decreases, and conversely.According to the above, and bearing in mind the spaces to which the unknowns belong (cf.( 20), (28), and ( 37)), we deduce that as the regularities of  and   increase, which means higher values for the exponents  and , the ones of the remaining unknowns decrease, that is  and  get smaller, and conversely.Consequently, no values yielding simultaneously either the least or the most regularity for each component of the solution are available, but only separately for each one of them.Certainly, the maximum or minimum regularity for a particular unknown in this latter case will not be achieved if the respective end of the corresponding interval is open.
We now go back to the well-definedness of s  by establishing the continuous inf-sup conditions for the bilinear forms   ,  P t1, 2u.While the corresponding proofs are similar to those of Lemma 2.7 from [26] and Lemma 3.6 of [11], and very close to that of Lemma 3.5 of [25], for sake of completeness we provide below the main details of them.
Proof.We begin by noticing that the values of  and  specified in (64) are compatible with the range r 2 `1 , 2 ´1 s required by Theorem 3.2 from [25], an existence result to be applied below.According to it, and since the pairs `1 ,  1 ˘and `2 ,  2 ˘result from each other exchanging  and , it suffices to prove (65) either for  " 1 or for  " 2. In what follows we consider  " 1, so that, given  P  1 :"   pΩq, we set   :" || ´2 , which belongs to   pΩq and satisfies ş Ω    " }} 0,;Ω }  } 0,;Ω .Thus, a straightforward application of the scalar version of Theorem 3.2 from [25] yields the existence of a unique  P  1, 0 pΩq such that ∆ "   in Ω,  " 0 on Γ.Moreover, the corresponding continuous dependence result reads }} 1,;Ω ď s   }  } 0,;Ω , where s   is a positive constant depending on .Next, defining  :" ∇ P L  pΩq, it follows that divpq "   in Ω, whence  P H  pdiv  ; Ωq ":  1 , and there holds }} 1 " }} ,div;Ω ď `1 `s   ˘}  } 0,;Ω .In this way, bounding by below with  :"  P  1 , and bearing in mind the definition of  1 (cf.(29b)) along with the foregoing identities and estimates, we arrive at (65) for  " 1 with  1 :" `1 `s   ˘´1 .The proof for  " 2 proceeds analogously, except for the fact that, given  P  2 :"   pΩq, and since  ă 2, one needs to define   :" As a consequence of Lemmas 4.3 and 4.4, and the boundedness properties given by ( 30)- (33), we are able to conclude now that the operator s  is well-defined.
Proof.Given  P  1 ˆ2 , a straightforward application of Theorem 2.1, Section 2.1 from [5] implies the existence of a unique p s , s q P  2 ˆ1 solution to (43).In turn, the a priori estimate provided in Corollary 2.1, Section 2.1, equation (2.15) of [5] establishes which, along with the aforementioned boundedness properties, yields (66) and ends the proof.
Similarly as for p  , and employing now ( [5], Cor.2.1, Sect.2.1, Eq. (2.16)), we observe that the a priori bound for the s  component of the unique solution to (43) reduces to which yields the same inequality as (66), but with a different constant, in particular depending additionally on s  1 .Therefore, as before, we still denote the largest of them by  s  , and simply say that the right hand-side of (66) constitutes the a priori estimate for both s  and s .

Well-definedness of the operator r 𝑇
In this section we employ the solvability result for perturbed saddle point formulations in Banach spaces provided by Theorem 3.4 of [15], along with the Banach-Nečas-Babuška Theorem (cf.[19], Thm.2.6), to show that, given an arbitrary p, vq P  2 ˆQ, equation ( 44) is well-posed for each  P t1, 2u, equivalently that   is well-defined.Since this result was already derived in Theorem 4.2 of [15] as an application of the abstract theory developed there, and more specifically of Theorem 3.4 from [15], we just discuss in what follows the main aspects of its proof.
To begin with, we introduce the bilinear forms ,  ,v : p  ˆ q ˆp  ˆ q Ñ  given by pp  ,   q, p  ,   qq :"   p  ,   q ` p  ,   q ` p  ,   q ´ p  ,   q, (68) and  ,v pp  ,   q, p  ,   qq :" pp  ,   q, p  ,   qq ´,v p  ,   q, for all p  ,   q, p  ,   q P   ˆ , and realize that (44) can be re-stated as: Find pĂ   , r   q P   ˆ such that  ,v ´´r   , r   ¯, p  ,   q ¯"   p  q ` p  q @p  ,   q P   ˆ .
In this way, the proof reduces to show first that the bilinear forms forming part of  satisfy the hypotheses of Theorem 3.4 from [15], and then to combine the consequence of this result with the effect of the extra term given by  ,v p¨, ¨q, to conclude that  ,v satisfies a global inf-sup condition.Indeed, it is clear from (38a), (38c), and the upper bound of   (cf.( 4)) that   and   are symmetric and positive semi-definite, which proves the assumption (i) of Theorem 3.4 from [15].Next, bearing in mind the definitions of   (cf.(38b)) and the spaces   and   (cf.(37)), and using again that   pΩq is isomorphic to the dual of   pΩq, we readily find that the null space   of the operator induced by   becomes and thus from which the assumption (ii) of Theorem 3.4 from [15], namely the continuous inf-sup condition for   , is clearly satisfied with constant r  :" s  ´1.In turn, while the continuous inf-sup condition for r   was already established in Lemma 2.9 from [26] (see also [15], Lem.4.1), for sake of clearness we provide below the main steps of its proof, which follows similarly to the one yielding the continuous inf-sup condition for b in the present Section 4.2.1.More precisely, given   P   :"   pΩq, we set  , :" |  | ´2   , which uses from (64) that  ě 2, and notice that there hold  , P   pΩq and ş Ω    , " }  } 0,;Ω } , } 0,;Ω .Then, we let   :" ∇ P L 2 pΩq, where  P  1 0 pΩq is the unique solution of the variational formulation: ş Ω ∇ ¨∇ " ´şΩ  ,  for all  P  1 0 pΩq, and deduce from the latter that divp  q "  , in Ω, which yields   P   :" Hpdiv  ; Ωq.In turn, denoting by  P the positive constant guaranteeing the Poincaré inequality:  P }} 2  1,Ω ď || for all p  ,   q P   ˆ , from which, under the assumption that, say Proof.Thanks to (77), (78), and the boundedness of   and   (cf.( 39), ( 40)), the unique solvability of (44) follows from a straightforward application of Theorem 2.6 from [19].In turn, the a priori estimate given by Theorem 2.6, equation (2.5) [19] reads which, along with the upper bounds for }  }  1  and }  }  1  derived from ( 39) and (40), yields (79) with  r  :" We end this section by observing from the definition of r  (cf.( 45)) and the priori estimates given by (79) for each  P t1, 2u, that for each p, vq P  2 ˆQ satisfying (76).

Solvability analysis of the fixed-point scheme
Knowing that the operators p  , s  , r  , and hence T as well, are well defined, we now address the solvability of the fixed-point equation ( 46).For this purpose, and in order to finally apply the Banach Theorem, we first derive sufficient conditions under which T maps a closed ball of  1 ˆ2 into itself.Thus, letting  be an arbitrary radius to be properly chosen later on, we define  pq :" t :" p 1 ,  2 q P  1 ˆ2 : }} 1ˆ2 ď u. (81) Then, given  P  pq, we have from the definition of T (cf.( 46)) and the a priori estimate for r  (cf.( 80)) that, under the assumption (cf.(76)) there holds For instance, defining letting  1 :" 2 0 , and imposing it is easily seen that (84) holds.We have therefore proved the following result.
Lemma 4.7.Assume that  and the data are sufficiently small so that there hold (84) and Then, T ` pq ˘Ď  pq.In particular, with the definition (85) of , and under the assumptions (86) and (87), the same conclusion is attained. We In turn, it is clear from (21d), and then subtracting and adding  to the factor  in the first term, that for each v P Q there holds pG , ´G, qpvq " ´1 tp 1 ´2 q p ´q `pp 1 ´1 q ´p 2 ´2 qq u ¨v, from which, proceeding as for the boundedness of G , (cf.( 23), ( 24)), that is employing the lower bound of  (cf.(4)), (9a), and the fact that } ¨}0,Ω ď |Ω| p´2q{2 } ¨}0,;Ω , we conclude that In this way, replacing (91) back into (90), we arrive at (88) and finish the proof.
which states that p s  ´s , s  ´s q P  2 ˆ1 is the unique solution of a problem like (43) with  " 0 and   ´ instead of   .In this way, proceeding as for the derivation of (66), which means applying the a priori estimate given by Corollary 2.1, Section 2.1, equation (2.15) from [5] (see also (67)), we find that Now, it is clear from (29d) that for each  P  2 there holds p  ´ qpq "  ´ pq " ´żΩ  tp 1 ´1 q ´p 2 ´2 qu, from which, applying Hölder's inequality, as we did for the boundedness of   (cf.( 30), ( 31)), and using that } ¨}0,;Ω ď |Ω| p´q{ } ¨}0,;Ω , we deduce that Finally, employing (95) in (94), we obtain (92) and conclude the proof.
It remains to prove the continuity of r  , which is provided by the following lemma.
Proof.Given p, vq and p, wq as indicated, we let, for each  P t1, 2u, r   p, vq :" r   P   and r   p, wq :" r   P   , where pr   , r   q P   ˆ and p r   , r   q P   ˆ are the corresponding unique solutions of (44), equivalently (cf.(70))  ,v ´´r   , r   ¯, p  ,   q ¯"   p  q ` p  q @p  ,   q P   ˆ , and  ,w ´´r   , r   ¯, p  ,   q ¯"   p  q ` p  q @p  ,   q P   ˆ .
It follows from ( 97) and (98), along with the definitions of the bilinear forms  ,v (cf.( 69)) and  ,v (cf.(38f)), that so that applying the global inf-sup condition (77) to pr   , r   q ´p r   , r   q, and then using (99) and the boundedness of  ,v (cf.( 39), ( 40)), we conclude that Having proved Lemmas 4.8-4.10,we now aim to derive the continuity property of the fixed point operator T. To this end, given ,  P  pq (cf.( 81)), we first recall from the definition of T (cf.( 46)) and Theorem 4.6 that, in order to define Tpq and Tpq, we need that the pairs `s  pq, p  p, s  pqq ˘and `s  pq, p  p, s  pqq satisfy (76).Then, according to the discussion at the beginning of the present section, we know that a sufficient condition for the latter is given by ( 84), which we assume in what follows.Alternatively, and as indicated there as well, (85) and ( 86) are in turn sufficient for (84).
Theorem 4.11.In addition to the hypotheses of Lemma 4.7, that is (84) and (87), or alternatively (85), (86), and (87), assume that Then, the operator T has a unique fixed point  P  pq.Equivalently, the coupled problem (41) has a unique solution p, uq P HˆQ, p, q P  2 ˆ1 , and p  ,   q P   ˆ ,  P t1, 2u, with  :" p 1 ,  2 q P  pq.Moreover, there hold the following a priori estimates Proof.We first recall that the assumptions of Lemma 4.7 guarantee that T maps  pq into itself.Then, bearing in mind the Lipschitz-continuity of T :  pq Ñ  pq (cf.( 104)) and the assumption (105), a straightforward application of the classical Banach theorem yields the existence of a unique fixed point  P  pq of this operator, and hence a unique solution of (41).Finally, it is easy to see that the a priori estimates provided by (52) (cf.Thm.4.1), (66) (cf.Thm.4.5), and (79) (cf.Thm.4.6) yield ( 106) and finish the proof.

The Galerkin scheme
We now introduce the Galerkin scheme of the fully mixed variational formulation ( 41), analyze its solvability by applying a discrete version of the fixed point approach adopted in Section 4.1, and derive the corresponding a priori error estimate.
Similarly to the remark right after (41) in Section 3.4, we highlight here that the second, fourth, and sixth rows of (107) constitute the discrete conservation of momentum properties, which are actually satisfied in an approximate sense.At the end of Section 6.2 we describe them explicitly in terms of suitable projection operators.

Discrete solvability analysis
In this section we proceed analogously to Sections 4.2 and 4.3 and establish the well-posedness of the discrete system (107) by means of the solvability study of the equivalent fixed point equation (113).In this regard, we emphasize in advance that, being the respective analysis very similar to that developed in the aforementioned sections, here we simply collect the main results and provide selected details of the corresponding proofs.
According to the above, we first aim to prove that the discrete operators p  ℎ , s  ℎ , and r  ,ℎ ,  P t1, 2u, and hence r  ℎ and T ℎ , are all well-defined, which reduces, equivalently, to show that the problems (108)-( 110) are well-posed.To this end, we now apply the discrete versions of Theorem 2.34 from [19], Theorem 2.1, Section 2.1 from [5], and Theorem 3.4 from [15], which are given by Proposition 2.42 of [19], Corollary 2.2, Section 2.2 from [5], and Theorem 3.5 from [15], respectively.More precisely, following similar approaches from related works (see, e.g., [11], Sect.4.2), our analysis throughout the rest of this section is based on suitable hypotheses that need to be satisfied by the finite element subspaces utilized in (107), which are split according to the requirements of the associated decoupled problems.Explicit examples of discrete spaces verifying these assumptions will be specified later on in Section 6.
We begin by addressing the well-definedness of p  ℎ , for which we let V ℎ be the discrete kernel of b, that is and assume that (H.1) there holds div `Hℎ ˘Ď Q ℎ , and (H.2) there exists a positive constant  d , independent of ℎ, such that sup Then, according to the definition of b (cf.(21b)), it follows from (114) and (H.1) that which says that V ℎ is contained in the continuous kernel V (cf.( 48)), and hence the discrete version of (50) is automatically satisfied, that is with  d "  :"   {.Recall here that   is the constant provided by inequality (49) with  " .In this way, it is clear from (117) that a satisfies the hypotheses given by Proposition 2.42, equation (2.35) from [19] with the constant  d , whereas (H.2) states that b fulfills ( [19], Prop.2.42, Eq. (2.36)) with the constant  d .We are thus in position to establish next the following result.
Similarly as observed for the continuous operator p  , we remark here that the right-hand side of (118) can also be assumed as the respective a priori estimate for p  ℎ .Furthermore, for the well-definedness of s  ℎ , we need to introduce the discrete kernels of  1 and  2 , namely and  2,ℎ :" t ℎ P  2,ℎ :  2 p ℎ ,  ℎ q " 0 @ ℎ P  2,ℎ u, respectively, and consider the following assumptions (H.3) there exists a positive constant s  d , independent of ℎ, such that sup sup (H.4) for each  P t1, 2u there exists a positive constant s  ,d , independent of ℎ, such that sup As a consequence of (H.3) and (H.4) we provide next the discrete version of Theorem 4.5.
Analogously as explained for the continuous operator s  , here we can also assume that, except for a constant  s  ,d depending additionally on s  1,d , the a priori estimate for s  ℎ , which follows now from Corollary 2.2, equation (2.25) of [5], is also given by the right-hand side of (123).
It remains to prove the well-definedness of r  ℎ :" p r  1,ℎ , r  2,ℎ q, for which we first observe that, being   and   symmetric and positive semi-definite in the whole spaces   and   , they certainly keep these properties in  ,ℎ and  ,ℎ , respectively, so that the assumption (i) of Theorem 3.5 from [15] is clearly satisfied.Next, given  P t1, 2u, we let  ,ℎ be the discrete kernel of   , that is  ,ℎ :" t ,ℎ P  ,ℎ :   p ,ℎ ,  ,ℎ q " 0 @ ,ℎ P  ,ℎ u, and consider the hypotheses (H.5) for each  P t1, 2u there holds div `,ℎ ˘Ď  ,ℎ , and (H.6) there exists a positive constant r  d ą 0, independent of ℎ, such that sup  ,ℎ P ,ℎ  ,ℎ "0 It follows from (124), the definition of   (cf.(38b)), and (H.5) that  ,ℎ :" t ,ℎ P  ,ℎ : divp ,ℎ q " 0u, whence, similarly to the case of p  ℎ ,  ,ℎ is contained in the continuous kernel   (cf.(71)) of   , thus yielding the discrete analogue of (72), that is In this way, it is clear from (127) that   satisfies the hypothesis (ii) of Theorem 3.5 from [15] with the constant r  d :" s  ´1, whereas (H.6) constitutes itself the corresponding assumption (iii).Consequently, a straightforward application of Theorem 3.5 from [15] implies the discrete global inf-sup condition for  (cf.(68)) with a positive constant r  ,d depending only on }  }, }  }, r  d , and r  d , and thus the same property is shared by   ℎ ,v ℎ for each p ℎ , v ℎ q P  2,ℎ ˆQℎ satisfying the discrete version of (76), that is We are now in position of establishing the well-definedness of r  ,ℎ for each  P t1, 2u.
Analogously to the continuous case, it follows from the definition of r  ℎ (cf.( 111)) and the a priori estimates given by (129) for each  P t1, 2u, that for each p ℎ , v ℎ q P  2,ℎ ˆQℎ satisfying (128).
Having established that the discrete operators p  ℎ , s  ℎ , r  ℎ , and hence T ℎ (under the constraint imposed by (128)), are all well defined, we now proceed as in Section 4.3 to address the solvability of the corresponding fixed-point equation (113).Then, letting  d be an arbitrary radius, we set  p d q :" t ℎ :" p 1,ℎ ,  2,ℎ q P  1,ℎ ˆ2,ℎ : and, reasoning analogously to the derivation of Lemma 4.7 (cf.beginning of Sect.4.3), we deduce that T ℎ maps  p d q into itself under the discrete versions of ( 84) and (87), which, denoting  0,d :" maxt1, for all  ℎ ,  ℎ P  p d q.
Consequently, we can establish next the main result of this section.
Theorem 5.4.Assume that  d and the data are sufficiently small so that (132) and (133) are satisfied, or alternatively that there holds (134), (135), and (133).Then, the operator T ℎ has a fixed point  ℎ P  p d q.

Preliminaries
In what follows we make use of the notations introduced at the beginning of Section 5.1.Thus, given an integer  ě 0, for each  P  ℎ we let P  pq and P  pq be the spaces of polynomials of degree ď  defined on  and its vector version, respectively.Similarly, letting x be a generic vector in   , RT  pq :" P  pq `P pqx and RT  pq stand for the local Raviart-Thomas space of order  defined on  and its associated tensor counterpart.In addition, we let P  p ℎ q, P  p ℎ q, RT  p ℎ q and RT  p ℎ q be the corresponding global versions of P  pq, P  pq, RT  pq and RT  pq, respectively, that is   p ℎ q :"  ℎ P  2 pΩq :  ℎ |  P   pq @ P  ℎ ( , P  p ℎ q :" v ℎ P L 2 pΩq : v ℎ |  P P  pq @ P  ℎ ( , RT  p ℎ q :" t ℎ P Hpdiv; Ωq :  ℎ |  P RT  pq @ P  ℎ u, and RT  p ℎ q :" t ℎ P Hpdiv; Ωq :  ℎ |  P RT  pq @ P  ℎ u.

Verification of the hypotheses (H.1)-(H.6)
We begin by observing from (155) that (H.1) is trivially satisfied, whereas (H.2) was proved in Lemma 5.5 from [14] (see, also, [8], Lem, 4.3) for the particular case given by  " 4 and  " 4{3.In turn, a vector version of (H.2) was established in Lemma 4.5 from [26] for  P p1, 2q in 2D (with local notation there given by  instead of ).In both cases, the preliminary result provided by Lemma 5.4 of [14] plays a key role in the respective proofs.While we could simply say, at least in 2D, that (H.2) follows basically from a direct extension of Lemma 4.5 from [26], we provide its explicit proof below for sake of completeness.To this end, following ( [26], Sect.4.1), we now introduce for each  P p1, `8q the space H  :"  P H  pdiv  ; Ωq Y Hpdiv  ; Ωq :  |  P W 1, pq @ P  ℎ ( , and let Π  ℎ : H  Ñ RT  p ℎ q be the global Raviart-Thomas interpolator (cf.[6], Sect.2.5).Then, we recall from Proposition 2.5.2 and equation (2.5.27) of [6] the commuting diagram property div `Π ℎ p q ˘"   ℎ `divp ˘˘@ P H  , where   ℎ :  1 pΩq Ñ   p ℎ q is the projector defined, for each  P  1 pΩq, as the unique element   ℎ pq P   p ℎ q such that ż In turn, it follows from Proposition 1.135 of [19] (see, also, [11], Eq. (A.5)) that there exists a positive constant   , independent of ℎ, such that for each  P p1, `8q there holds On the other hand, while here we could use again ( [14], Lem.5.4), we prefer to resort to the slightly more general result provided by Lemma A.2 of [11], thus giving a greater visibility to it, which establishes that, given an integer  such that 1 ď  ď  `1, and given ,  P p1, `8q, such that  ď  ď  ´ if  ă , or  ď  ă `8 if  " , there exists a positive constant , independent of ℎ, such that Note that for the first set of constraints on  and , there holds   ´  ě ´1, which yields  `  ´  ě 0, whereas for the second one, there holds  `  ´  "  ´1 `  ě   , thus proving that in any case the power of ℎ in (159) is non-negative.In this way, it follows from (159) that, for  " 1, and under the specified ranges of  and , there exists a positive constant  Π , independent of ℎ, such that (cf.[11], Lem.A.3) In particular, for  ă  and  " 2, the inequality  ď  ´ becomes  ě 2 `2 , so that for the resulting range of , that is  P " 2 `2 , 2 ˘in 2D, and  P Analogue identities and inequalities to those stated above are valid with the tensor and vector versions of Π  ℎ and   ℎ , which are denoted by Π  ℎ and   ℎ , respectively.We are now in position to prove that (H.2) holds.Lemma 6.1.Under the ranges for  and  specified by (64), there exists a positive constant  d , independent of ℎ, such that Proof.Given v ℎ P Q ℎ , v ℎ " 0, we set v ℎ, :" |v ℎ | ´2 v ℎ , which belongs to L  pΩq, and notice that ż Ω v ℎ ¨vℎ, " }v ℎ } 0,;Ω }v ℎ, } 0,;Ω .
Next, we let  be a bounded convex polygonal domain that contains s Ω, and define g :" It is readily seen that g P L  pq and }g} 0,; " }v ℎ, } 0,;Ω .Then, applying the elliptic regularity result provided by Corollary 1 from [21], we deduce that there exists a unique z P W 2, pq X W 1, 0 pq such that: ∆z " g in , z " 0 on B.Moreover, there exists a positive constant  reg , depending only on , such that }z} 2,; ď  reg }g} 0,; "  reg }v ℎ, } 0,;Ω .
In turn, noting that the range for  (cf.( 64)) fits into the one for  in (161), we can apply this inequality (with  " ) and the regularity estimate (165), to arrive at which, combined with (167), implies Consequently, bounding below the supremum in (162) with  ℎ , and making use of (166), the analogue of ( 157), (163), and (169), we conclude the required discrete inf-sup condition with the constant  d :" ` `Π  reg ˘´1 .
Furthermore, for the hypotheses (H.3) and (H.4), we first stress that (H.3) corresponds exactly to (H.5) of [11], and hence we omit most details and refer to Section 5.2, Lemma 5.2 from [11].We just make a few remarks here.First of all, we observe that the discrete kernels of the bilinear forms  1 and  2 coincide algebraically, which reduces to   ℎ :" t ℎ P RT  p ℎ q : divp ℎ q " 0 in Ωu.Then, we let Θ  ℎ : L 1 pΩq Ñ   ℎ be the projector defined similarly to (157), that is, given In this way, a quasi-uniform boundedness property of Θ  ℎ in 2D (cf.[11], Eq. (5.8)), along with the properties of the operators   (cf.Lem.4.2), play a key role in the proof of (H.3).Whether the aforementioned boundedness is satisfied or not in 3D is still an open problem, and hence, similarly to [11], the assumption (H.3) is the only aspect of the analysis in this section that does not hold in 3D.All the other conditions are valid in both 2D and 3D.Regarding (H.4), we remark that the discrete inf-sup conditions for  1 and  2 , which adapt the continuous analysis from Lemma 4.4 to the present discrete setting, follow from slight modifications of the proofs of Lemma 4.5 from [26] and Lemma 5.3 from [11].Further details are omitted here.
Finally, it is clear from (155) that (H.5) is trivially satisfied, whereas (H.6) was proved precisely by Lemma 4.5 from [26].Alternatively, for the discrete inf-sup condition for   we can proceed analogously to the proof of Lemma 6.1 by observing that the range of  (cf.(64), recall that   :" Hpdiv  ; Ωq) also fits into the one for  in (161), whence this inequality can be applied to  "  as well.
From the last columns we see that the potential and transport balance equations are satisfied to machine precision while the error for the momentum balance is higher.This may be explained by the presence of the term  ℎ on the right-hand side (which has a H(div)-component).
Example 2. In addition, and in order to illustrate the implementation of fixed-point solvers, we have realized numerically Picard versions of the linearization of (107).In case A we follow the fixed-point structure used in the analysis of Section 5.1, that is, solving sequentially problems (108) Ñ (109) Ñ (110), and iterating until the ℓ 2 -norm of the vector containing the residual of the Picard iterates reaches 10 ´8.Next, in case B we choose a different fixed-point splitting where we apply two modifications with respect to case A. First, in (110) instead of the linear functional for the second discrete electrostatic potential equation (discrete version of (29d)) we consider p ℎ q :" ´şΩ   ℎ and the coupling term appears as a bilinear form contribution (and no longer as part of the linear functional), say p p ℎ , p 1,ℎ ,  2,ℎ qq :" ż Ω  ℎ p 1,ℎ ´2,ℎ q.
We focus on the number of Picard iterations required in each case, displaying the obtained results in Table 2.
While we confirm that all methods give exactly the same errors (and consequently also the same convergence rates, which are optimal in view of the theoretical bounds), from the number of fixed-point iterations we readily note that case B performs much better than case A, for the two polynomial degrees we tested  " 0,  " 1.This behaviour could be explained by the stability of different linearizations of advective nonlinearities and by the strength of the coupling for this particular choice of model parameters.We stress that the analysis of case B is, however, not at all straightforward since the decoupled linear electrostatic potential problem resulting from the first modification is no longer symmetric.For sake of reference we also tabulate total errors and number of nonlinear iterates obtained with the method we use also in Examples 1 and 3: an exact Newton-Raphson linearization (labelled here as case C).Needless to say, the latter is actually the one that one could eventually employ in practical computations.Samples of the approximate solutions (only the mixed variables) computed with the method in case A are portrayed in Figure 3.
Example 3. We conclude this section with an application problem where we demonstrate the use of the mixed finite element scheme in simulating the transport process in an electrokinetic system with an ion-selective interface, where the development of an electroosmotic instability is expected.The problem configuration is adopted from [17,18].This system corresponds to a transient counterpart of (8) in the absence of external forces and sources (f " 0,  "   " 0), where the following additional terms appear in the momentum and concentration equations (note also the different scaling of  on the right-hand side of the momentum balance, required to match the adimensionalization in [18]) ´1 Sc B  u ´divpq " p 1 ´2 q 1 2 2 , ´B   `divp  q " 0.
The time derivatives are discretized using backward Euler's method.In the problem setup a boundary layer is present in the vicinity of the solid boundary (the bottom edge of the rectangular domain), and therefore we employ a graded mesh with a higher refinement close to the layer.For this problem we select the second-order family of finite element subspaces (setting  " 1 in Sect.6.1), which gives for the chosen mesh 865 201 degrees of freedom.
The physical properties of the system are as follows.The cation species is Na `having diffusivity  1 " 1 and the anion species is Cl ´with the same diffusivity  2 " 1.The dynamic viscosity of the mixture is  " 1.Initial conditions are given by u " 0, and a 2% random perturbation on a linearly varying initial ionic concentrations  1 " p2 ´q,  2 " , where  is a uniform random variable between 0.98 and 1.On the top boundary we set  1 "  2 " 1, u " 0, and an applied voltage of  " 120.On the bottom boundary we impose  " 0,  1 " 2,  2 ¨ " 0, and u " 0. On the vertical walls we prescribe periodic boundary conditions.The other model parameters take the values  " 8 ˆ10 ´6, Sc =10 3 , and we use a timestep ∆ " 10 ´6.We plot snapshots of

Lemma 4 . 4 .
For each  P t1, 2u there exists a positive constant s   such that sup P   "0

Figure 2 .
Figure 2. Example 1. Error history associated with the finite element family (155) with  " 0 in 3D for primal variables (top left) and mixed variables (top right, including the velocity gradient using the postprocess in (153)), and samples of approximate primal variables (velocity streamlines u ℎ , iso-surfaces of postprocessed pressure  ℎ , electrostatic potential  ℎ , and positive ion concentration  1,ℎ ; bottom plots).In all mesh refinements the number of Newton-Raphson iterations was 4.

Figure 3 .
Figure 3. Example 2. Samples of approximate mixed variables (stress magnitude, electric field magnitude and arrows, and ionic fluxes) obtained with the fixed-point algorithm labelled case A, and for  " 1.

Figure 4 .
Figure 4. Example 3. Samples of approximate velocity (top) and anion concentration (bottom) at times  " 10 ´4 and 10 ´3 (left and right, respectively), produced with the mixed method and using  " 1.

Table 2 .
Example 2.Total error, experimental rates of convergence, and number of iterations required for two types of fixed-point methods as well as for Newton-Raphson linearization.