On a vorticity-based formulation for reaction-diffusion-Brinkman systems

We are interested in modelling the interaction of bacteria and certain nutrient concentration within a porous medium admitting viscous flow. The governing equations in primal-mixed form consist of an advection-reaction-diffusion system representing the bacteria-chemical mass exchange, coupled to the Brinkman problem written in terms of fluid vorticity, velocity and pressure, and describing the flow patterns driven by an external source depending on the local distribution of the chemical species. A priori stability bounds are derived for the uncoupled problems, and the solvability of the full system is analysed using a fixed-point approach. We introduce a primal-mixed finite element method to numerically solve the model equations, employing a primal scheme with piecewise linear approximation of the reaction-diffusion unknowns, while the discrete flow problem uses a mixed approach based on Raviart-Thomas elements for velocity, Nedelec elements for vorticity, and piecewise constant pressure approximations. In particular, this choice produces exactly divergence-free velocity approximations. We establish existence of discrete solutions and show their convergence to the weak solution of the continuous coupled problem. Finally, we report several numerical experiments illustrating the behaviour of the proposed scheme.


1.
Introduction. Reaction-diffusion systems can explain many phenomena taking place in diverse disciplines such as industrial and environmental processes, biomedical applications, or population dynamics. For instance, non-equilibrium effects associated to mass exchange and local configuration modifications in the concentration of species are usually represented in terms of classical reaction-diffusion equations. These models allow to reproduce chaos, spatiotemporal patterns, rhythmic and oscillatory scenarios, and many other mechanisms. Nevertheless, in most of these applications the reactions do not occur in complete isolation. The species are rather immersed in a fluid, or they move within (and interact with) a fluid-solid continuum, and the motion of the fluid inevitably affects that of the species. In some circumstances, reciprocal effects might be substantially large, therefore leading to local changes in the observed flow patterns. More specifically, here it is assumed that the medium where the chemical reactions develop is a porous material saturated with an incompressible fluid. We remark that we are interested in viscous flows and consequently, the linear momentum and mass conservation for the fluid are governed by the Brinkman equations. The regime under study is relevant to fluids within highly porous structures composed by light fixed particles [7]. As indicated above, we also suppose that the local fluctuations of a species' concentration is important enough to affect the fluid flow. In turn, the reaction-diffusion equations include additional terms accounting for the advection of each species with the fluid velocity, therefore improving the mixing and interaction properties with respect to those observed under pure diffusion effects [13]. Diverse types of continuum-based models resulting in reaction-diffusion-momentum equations have been applied for e.g. enzyme reactions advected with a known Poiseuille velocity profile [16], population kinetics on moving domains [28,35], the control of drug release within soft living tissue [9,15], or biochemical interactions on growing surfaces [12].
We aim at developing numerical solutions to a class of similar systems, also addressing the solvability of the governing equations, the expected properties of the underlying solutions, and the convergence behaviour of suitable finite element schemes. More precisely, at both continuous and discrete levels, the Brinkman equations are set in a mixed form (that is, the associated formulation possesses a saddle-point structure involving the vorticity as additional unknown) whereas the formulation of the reaction-diffusion system is written exclusively in terms of the primal variables, in this case the species' concentration. Such a structure is motivated by the need of accurate vorticity rendering, obtained directly without postprocessing it from typically low-order discrete velocities (which usually leads to insufficiently reliable approximations). The nonlinear reaction-diffusion system will be placed in the framework of general semilinear parabolic equations and its wellposedness, for a fixed velocity, follows from the classical theory of Ladyženskaja [22]. On the other hand, for a fixed species' concentration, a number of recent Brinkman formulations based on vorticity, pressure and velocity are available from the literature [3,5,6,19] (see also [14]). Here we adopt the one introduced in [4], where the well-posedness of the flow equations is established thanks to the classical Babuška-Brezzi framework. In turn, the fully coupled problem is analysed using a Schauder fixed-point approach. Perhaps the closest contributions in the spirit of the present study are the error estimation for finite element formulations to doubly-diffusive Boussinesq flows presented in [2], the two-sidedly degenerate chemotaxis-fluid coupled system from [11] whose solvability relies on a regularisation step combined with compactness arguments; and the general analysis for degenerate parabolic equations idealising reactive solute transport in porous media recently carried out in [21].
A further goal is to propose a numerical method based on piecewise linear and continuous polynomials for the species' concentrations, whereas the set of flow equations is discretised using a mixed approach based on Raviart-Thomas elements for the approximation of the velocity, Nédélec elements for vorticity, and piecewise constants for the pressure. The computational burden of this algorithm is comparable to classical low-cost approximations as the so-called MINI element for velocitypressure formulations. On top of that, the proposed finite element method provides divergence-free velocities, thus preserving an essential constraint of the underlying physical system. Based on the present formulation, the stability properties of different staggered coupling techniques have been recently analysed in [24]. Other recent contributions regarding the numerical discretisation for related models for reactive solute transport in porous media (but using Darcy descriptions of flow), include the formulation and convergence studies for mixed and conformal methods as presented in [1,8,20,27,30,31].
After this introductory section, the remainder of the paper is structured in the following manner. Section 2 states the model problem and introduces the concept of weak solutions, together with an auxiliary coupled problem. The solvability of these auxiliary equations is established in Section 3, using a fixed-point approach. The proof of uniqueness of weak solution is postponed to Section 4. Section 5 deals with the numerical approximation of the model problem, using mixed finite elements and a first-order backward Euler time advancing scheme. We provide in Section 6 a set of numerical tests to illustrate the properties of the numerical scheme and the features of the coupled model, and we close in Section 7 with a summary and discussion of possible extensions to this work.

Problem statement.
Let Ω ⊂ R 3 , be a simply connected, and bounded porous domain saturated with a Newtonian incompressible fluid, where also a bacterium and some chemical species (diffusible agents or nutrients) are present. Viscous flow in the porous medium is usually modelled with Brinkman equations stating momentum and mass conservation of the fluid. In addition, the Reynolds transport principle applied to mass conservation of the interacting species yields an advection-reaction-diffusion system. The physical scenario of interest can be therefore described by a coupled system written in terms of the fluid velocity u, the rescaled fluid vorticity ω, the fluid pressure p, and the volumetric fraction (or concentration) of the bacterium c and of the chemical substance s. Let t denote the time variable taking values in the interval (0, T ], where T is a given final time. Then for a.e. (x, t) ∈ Ω T := Ω × (0, T ], we consider where µ is the fluid viscosity (in the considered regime, it is assumed independent of the bacteria and chemical concentrations), K(x) is the permeability tensor rescaled with viscosity, sg represents the force exerted by the bacteria on the fluid motion (where g is the gravitational acceleration, or any other particular force to which the species comply), and f (x, t) is an external force (e.g. centrifugal) applied to the porous medium. Moreover, D c , D s , G c , G s are concentration-dependent coefficients determining, respectively, the nonlinear diffusivities and reaction kinetics (representing production and degradation) of bacteria and chemical species. Equations (2.1) are complemented with the following boundary and initial data: representing that the species cannot leave the medium and that a normal velocity u ∂ together with a compatible vorticity trace ω ∂ are imposed along the domain boundary ∂Ω. The presentation will be restricted however to homogeneous boundary conditions. The analysis readily extends to the non-homogeneous case by classical lifting arguments.
Our first main result states existence of weak solutions to the reaction-diffusion-Brinkman equations. The analysis of uniqueness will be postponed to Section 4. The proof of this result is based on an application of Schauder's fixed-point theorem (in an appropriate functional setting) with the derivation of a priori estimates and the compactness arguments, to the following uncoupled system ∂ t c + u · ∇c − div(D c (c)∇c) = G c (c, s), ∂ t s + u · ∇s − div(D s (s)∇s) = G s (c, s), whereŝ is a fixed function.
3. Existence result: Fixed-point approach. In this section we prove the existence of solutions to (2.7), looking first at the solvability of the uncoupled systems. Theorem 3.1. Let (X , ·, · X ) be a Hilbert space. Let A : X × X → R be a bounded symmetric bilinear form, and let G : X → R be a bounded functional. Assume that there existsβ > 0 such that Then, there exists a unique x ∈ X , such that Moreover, there exists C > 0, independent of the solution, such that Let us also consider the kernel of the bilinear form Ω q div u dx, that is Moreover, we endow the space H(curl; Ω) with the following µ-dependent norm: ,Ω , and recall the following inf-sup condition (cf. [17]): there exists β 2 > 0 only depending on Ω, such that The fixed-point method. In view of invoking Schauder's fixed-point theorem to establish solvability of (2.7), we introduce the following closed subset of the Banach space L 2 (Ω T ): and Lemma 3.4, respectively. Our idea is quite simple, for a givenŝ ∈ K s we find first the velocity u and then the pair (c, s) solutions of (2.7). Next, and thanks to Theorem 3.1, we can assert the solvability of the Brinkman equations for a fixed s ∈ K s and for any t > 0.
As the next step, we consider a constant λ > 0 and define the auxiliary variables (φ c , φ s ) by setting c = e λt φ c and s = e λt φ s . We now introduce a map R : K s → K s such that R(ŝ) = s, where s solves (3.5). The goal is to prove that such map has a fixed-point. First, let us show that R is a continuous mapping. Let (ŝ n ) n be a sequence in K s andŝ ∈ K s be such that s n →ŝ in L 2 (Ω T ) as n → ∞. Let us then define s n = R(ŝ n ), i.e., s n is the solution of (3.5) associated withŝ n and the solution u of (3.3). The objective is to show that s n converges to R(ŝ) in L 2 (Ω T ). We start with the following lemma: Lemma 3.4. Let (c n , s n ) n be the solution to problem (3.5) and recall thatŝ ∈ K s . Then Proof. (i) We replace the strong form of (3.5) by Multiplying (3.7) by −c − n = c n − |c n | 2 and integrating over Ω, we obtain Now, we use that div u = 0 and (2.5) to deduce According to the positivity of the second and the fourth terms in the left-hand side of (3.8), we get 1 2 Since c 0 is nonnegative, we deduce that c − n = 0, and reasoning similarly we have that s − n = 0. Next, we let k c , k s ∈ R such that k c ≥ c 0 L ∞ (Ω) and k s ≥ s 0 L ∞ (Ω) . Let us consider t ∈ (0, T ), and proceed to multiply (3.7) by ξ c := (c n − e (β−λ)t k) + , with β ≥ λ, and with k = sup{k c , k s }, and to integrate over Ω, which yields that there exists some constant C 6 > 0 such that where ξ s := (s n − e (β−λ)t k) + . Following an analogous treatment for s n , we get for some constant C 7 > 0. We next observe that which inserted into (3.9) and (3.10), implies that Finally for β ≥ λ ≥ C 6 + C 7 , an application of (3.11) yields and exploiting that c 0 , s 0 ≤ k in Ω, we conclude that c n (·, t) ≤ e (β−λ)t k and s n (·, t) ≤ e (β−λ)t k in Ω × (0, T ).
(ii) We multiply equation (3.7) by c n and integrate over Ω to obtain ≤ Ω e −λt G c (e λt c n , e λt s n )c n dx. (3.12) Invoking the boundedness of c n and s n , we get that the integral on the right-hand side is bounded independently of n. Then using that we deduce from (3.12) the following bound valid for some constant C 8 > 0 independent of n. This completes the proof of (ii).
(iii) Next, we multiply (3.7) by ϕ ∈ L 2 (0, T ; H 1 (Ω)) and use the boundedness of c n and s n to get for some constants C 9 , C 10 , C 11 > 0 independent of n. Consequently, we end up with the bound Working on the same lines as for (c n ) n , the counterpart of (ii) − (iii) and (3.13) for (s n ) n can be obtained. Finally, (iv) is deduced from (ii) and (iii).

4.
Uniqueness of weak solutions. The following result completes the analysis of unique solvability for the reaction-diffusion-Brinkman system. , f = f 2 , respectively. Then, for any t ∈ (0, T ] there exists C > 0 such that In particular, there exists at most one weak solution to the reaction-diffusion-Brinkman system (2.1)-(2.2).

Numerical approximation.
In this section, we construct a primal-mixed fullydiscrete scheme for the approximation of the coupled system (2.1). The finite element spaces characterising the semi-discrete problem are then specified, and a proof of convergence to the unique weak solution (presented in Definition 2.1) is outlined.
5.1. The semi-discrete scheme. Let T h be a regular family of triangulations of Ω by tetrahedra K of maximum diameter h. Given an integer k ≥ 0 and S ⊂ R 3 , by P k (S) we denote the space of polynomial functions defined in S of total degree up to k, and define the following finite element subspaces where for any K ∈ T h (Ω), the local lowest order Raviart-Thomas and edge space of Nédélec type, are defined as RT 0 (K) = P 0 (K) 3 ⊕ P 0 (K)x, and ND 1 (K) = P 0 (K) 3 ⊕ P 0 (K) 3 × x, respectively. Then, a Galerkin semi-discretisation associated to the formulation introduced in Definition 2.1 reads:

Euler time discretisation.
Let c 0 h = P V h (c 0 ), s 0 h = P V h (s 0 ) be appropriate projections of the initial data, and consider the following fully discrete method arising after backward Euler time discretisation using a fixed time step ∆t = T /N : where the forcing term in the momentum equation is discretised explicitly. We also stress that, owing to the choice (5.1), the discrete velocities u h generated using either (5.2) or (5.3) are exactly divergence-free, that is, div u h (t) = 0 and div u n h = 0 in Ω (cf. [4]).

5.3.
Sketched convergence analysis. Note that the positive and negative cuts test functions may not be admissible in finite element context. Therefore, we cannot obtain the maximum principle as in the continuous case. Comparing to the continuous case, here we assume that the initial condition (c 0 , s 0 ) is only in L 2 (Ω, R 2 ). In order to derive stability estimates, we employ the same kind of test functions in the same combination of equations used in Section 3; whereas the chain rule for time derivatives is now replaced by the convexity inequality a n (a n − a n−1 ) ≥ |a n | 2 2 − a n−1 2 2 , to treat the finite difference discretisation of the time derivatives. Proceeding in this way, we can obtain estimates (uniform in both h and ∆t) for the discrete Brinkman solutions u h L 2 (0,T ;H(div;Ω)) + ω h L 2 (0,T ;H(curl;Ω)) + p h L 2 (0,T ;L 2 0 (Ω)) ≤ C, (5.4) and for the advection-reaction-diffusion system for some constant C > 0. The next goal is to establish the relative compactness in L 2 (Ω T ) of the sequences (c h , s h ), which is achieved by constructing space and time translates and using the a priori estimates given above.
for all r ∈ R 3 and for all τ Proof. Let us introduce the space translates for all m c h , m s h ∈ V h . We integrate (5.9) with respect to time σ ∈ [t, t + τ ] (with 0 < τ < T ) and consider m c h = T τ c h and m s h = T τ s h as test functions in the resulting equations. Therefore Now, we examine these integrals separately. For I 1 we have and similarly |I 2 | ≤ C τ , for some constant C > 0. Herein we have used the Fubini theorem (recall that t+τ t dσ = τ = σ σ−τ dt), the divergence-free condition for the discrete velocity, the Hölder inequality and the L 2 −bounds for (c h , s h ), (∇c h , ∇s h ) and (∇T h c h , ∇T h s h ). Keeping in mind the growth assumptions on G c , G s , we can then apply Cauchy-Schwarz inequality to deduce that |I 3 | + |I 4 | ≤ C τ, for some constant C > 0. Collecting these inequalities we readily get Furthermore, the definition of (c h ,s h ) together with (5.5) eventually implies that which establishes (5.7), therefore finishing the proof.
Note that as a consequence of (5.4)-(5.5), Lemma 5.1, and Kolmogorov's compactness criterion (cf. [10, Theorem IV.25]), we can assert that there exists a subsequence of (c h , s h , u h , ω h , p h ), not re-labeled, such that, as h → 0, (c h , s h ) → (c, s) strongly in L 2 (Ω T , R 2 ) and in L 2 (0, T ; H 1 (Ω, R 2 )) weakly, in L 2 (0, T ; H(div; Ω)) × L 2 (0, T ; H(curl; Ω)) × L 2 (0, T ; L 2 0 (Ω)) weakly. These convergences allow us to identify the limit (c, s, u, ω, p) as the weak solution of (2.1), and the convergence result is summarised as follows. We observe that the error generated by the fully discrete scheme (5.3) has two components: one due to the spatial discretisation and depending on the meshsize h, and the error due to the time discretisation depending on the timestep ∆t. Given the approximation properties of the employed finite element spaces and the time stepping method (see e.g. [37]), we can expect the following convergence rates for the proposed method (c(·, t n ), s(·, t n ), u(·, t n ), ω(·, t n ), p(·, t n )) − (c n h , s n h , u n h , ω n h , p n h ) ≤ C 1 h + C 2 ∆t, (5.10) with t n = n∆t, for n = 1, . . . , N and C 1 , C 2 > 0 are constants independent of h and ∆t. Here (c n h , s n h , u n h , ω n h , p n h ) denotes the sequence generated by (5.3) for all n = 1, . . . , N . A rigorous derivation of this space-time error estimate is part of ongoing developments.
6. Numerical tests. We finally present a set of elementary examples to illustrate the properties of the model and of the proposed finite element method. The coupling between Brinkman and reaction-diffusion equations will be implemented either via: a) a fully monolithic solution, b) an iterative Picard method splitting linear Brinkman and fully explicit reaction-diffusion equations, and c) an iterative Picard method with an embedded Newton algorithm for the linearisation of the reactiondiffusion system. Rigorous convergence proofs for the Newton iteration and various linear iterative schemes applied to similar problems in non-viscous flow in porous media can be found in e.g. [27,34]. For our particular scenario, efficient coupling   Figure 1. Example 1. Convergence tests for the spatial (left) and temporal (right) discretisation via mixed P 1 × P 1 × RT 0 × P 1 × P 0 finite elements and backward Euler time stepping applied to (2.1).
strategies and thorough comparisons between dedicated solvers have been reported in [24].
The model parameters are set as µ = 1, D c = 0.01, D s = 0.5, g = (0, −1) T , and K = I; the reaction kinetics are specified by whereas f is computed using (6.1) and the momentum equation in (2.1). Note that these solutions satisfy the mass conservation equation, while boundary and initial conditions (2.2) are imposed according to (6.1). Also, a non-homogeneous source term is incorporated on the right hand side of the reaction-diffusion equations and constructed using the manufactured solutions. Errors associated to splitting algorithms are avoided by invoking the exact full Jacobian and performing Newton iterations until convergence at each time step.
The mesh convergence is investigated by fixing ∆t =1.6e-4 and computing errors between the exact solutions from (6.1) and approximations obtained with our mixedprimal method on a sequence of successively refined structured triangulations of Ω, at the final adimensional time T =5e-3. That is, for each field η, we compute e h (η) := η N − η N h , where η N = η(·, T ), and · is either the H 1 , H(div), or L 2 −norm. On the other hand, the error history associated to the time discretisation is studied by fixing a small meshsize h = 2.43e-5 and computing accumulative errors up to T = 0.1, and decreasing ∆t. That is, we measure errors in the ∞ (0, t; ·)−norm: E ∆t (η) := N n=0 η n − η n h (·) . Figure 1 depicts a convergence analysis for the propose method with respect to meshsize and timestep. All plots indicate an overall first order convergence in both space and time, as expected from the involved approximations yielding (5.10). In Table 6 we provide the iteration count to reach convergence of the Newton method for the first test of spatial accuracy. The values represent the average of required steps at each resolution, and indicate that the convergence is independent of the refinement level.
Example 2. Bacterial bioconvection. Next we focus on the interaction between bacteria (with concentration c) and oxygen (described by s) in a small chamber, similar to the case studied in e.g. [23]. Let us define the function r(s) , where s * is an oxygen threshold, below which the chemotactic convection is turned off. We consider a semi implicit approximation of the cross diffusion, the reaction part for the oxygen conservation equation is given by G s (c, s) = βcr(s), and the remaining model functions and adimensional parameters are set as follows: K = K 0 + K 0 η(x) (with K 0 = 7700 and η a uniformly distributed random field), g = (0, −γ) T , s * = 0.3, D c = 0.01, D s = 0.25, ε = h. The computational domain is a disk of radius 1 2 centred at ( 1 2 , 1 2 ), discretised into an unstructured mesh of 13972 points and 27942 triangular elements, and a timestep of ∆t =1e-3 is employed. For this example we solve via fixed point iterations the coupling between the Brinkman problem and the set of reaction-diffusion equations. The system is initially at rest (zero velocity, vorticity, and pressure), having a concentration of bacteria near the top of the disk c 0 = 1−(1+exp(−50 (x − 0.5) 2 + (y − 0.9) 2 )) −1 , and an homogeneous oxygen content s 0 = 1. In Figure 2 we show three snapshots (at advanced time) of the obtained numerical solutions on different regimes characterised by α, β, γ. In all cases we observe the species heading to the bottom of the disk, generating vortical flow around the zones of high concentrations of oxygen and bacteria. The plots indicate that transitional flow occurs (mainly due to the triggering of unstable modes), which can be captured in a robust manner by the numerical method, even in the most unstable regime.  set a timestep of ∆t =2.5e-3, and run the simulation until T = 6. First, we use the method of characteristics to treat the advective terms, whereas the reaction kinetics are considered explicit. The obtained numerical solutions are depicted in Figure 3. We investigate the convergence properties of the coupling schemes. In a second round of tests we now use a Picard iterative iterative scheme to couple the Brinkman and the FitzHugh-Nagumo equations, and consider the reaction terms implicitly, so that an embedded Newton iteration is applied inside each fixed-point solution. Figure 4 portrays the evolution of the number of required outer fixed-point and inner Newton steps to achieve convergence (set with a tolerance of tol=1e − 6 for both schemes). We see that the number of Picard steps needed to converge are higher at the beginning of the simulation and they slowly decrease, whereas the Newton steps applied to the nonlinear reaction terms do not change substantially. We also perform the 3D counterpart of this example, defined in the polygonal domain Ω = (0, 1.8) × (0, 1) × (0, 0.6) filled with 60 large particles (of radii 1.5e-3 and randomly distributed on Ω) with a permeability 100 times higher than in the rest of the porous matrix. The forcing term is specified by g = d 2 (1, 1, 1) T and f = 0, whereas the conductivity matrix and remaining parameters are chosen as above. The domain is discretised into an unstructured mesh of 126034 tetrahedral elements sharing 21952 nodes, and the same timestep as above is employed. All sides of the box are provided with zero-fluxes for voltage and recovery fields, zero tangential vorticity, and zero normal velocities, and as initial condition we excite the bottom left part of the domain. For this test a Newton algorithm is used to linearise the reaction-diffusion equations, and outer Picard iterations are applied to couple that system together with the set of Brinkman equations. The obtained numerical solutions are depicted in Figure 6. As in its 2D counterpart, in this test we notice that a propagating front for the potential moves towards the positive x−axis, followed by the slower recovery variable front. Due to the heterogeneity of    conductivities and permeabilities, preferential velocity patterns start to form. We also notice that vorticity (in this case, we show only its magnitude) clearly marks the regions of contact between high gradients of potential and recovery field.  Example 4. Intracellular calcium-induced calcium release. In closing this section we present a simulation of the interaction between two species c, s (the concentrations of cytosolic and sarcoplasmic calcium, respectively) inside a cardiac cell. This phenomenon has been studied in terms of the reacting species alone (cf. [18,26]), and also including the active cell contraction while solving for the underlying finite elasticity equations [29,38]. In contrast, here we assume that the species interact with an interstitial fluid occupying the sarcoplasmic reticulum, which is in turn pictured as a porous medium with a non-homogeneous permeability distribution. A 2D geometry reconstructed from confocal images is used and a triangular mesh of 46610 elements and 23741 vertices is generated. The time advancing algorithm uses a fixed time step of ∆t = 0.01. The process consists in opening 80 channels of cytosolic calcium located randomly within the myocyte, and observing how these propagate thorough plasma into the whole cell. We also consider that the permeability is higher in the vicinity of these channels. The minimal reaction-diffusion model proposed in [18] involves the following specialisation for constant and speciesdependent coefficients: where η is a step function assuming the value 1 on disks centred at uniformly distributed random locations corresponding to the channels, and zero elsewhere. The influence of calcium patterns into the flow behaviour of the plasma is encoded in the forcing term for the momentum equation f (s) = γ 0 |s|(f 0 ⊗ f 0 )∇s, which might be regarded as the flow-counterpart to the active force proposed in [28].
Starting from the open channels, the cytosolic calcium starts to propagate towards other channels, also following a higher diffusion in the preferential direction f 0 (aligned with the x−axis).
7. Concluding remarks. The interaction of viscous flow in porous media and reaction-diffusion phenomena occur in a variety of applicative scenarios. Here we have advocated the solvability analysis and the numerical approximation of a special class of problems related to Brinkman and nonlinear advection-diffusion-reaction equations. The analysis of the coupled set of equations is framed under the theory of fixed-point schemes combined with an abstract result for saddle-point problems.
A mixed-primal finite element method has been proposed and we have discussed its stability and convergence properties. We have provided numerical validation of the spatio-temporal convergence of the method, and have also addressed the simulation of a number of biochemically-oriented examples. Our theoretical framework currently does not account for higher order terms nor the chemotactical convection simulated in Example 2, and further investigation is needed in this direction. We also need to cover the dependence of the momentum force on the gradient of the chemical concentrations (studied numerically in Example 4). Regarding other model generalisations, interesting topics be explored include cross-diffusion effects, nonlinear convective fluxes, a non-laminar flow regime, viscosity depending on the chemical concentrations, or variable density flows. Finally, we mention that the study of mixed formulations involving embedded saddle-point problems and the convergence analysis of the associated fully-mixed finite element schemes can be the subject of suitable extensions to the mathematical aspects of this contribution.