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

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 well-posedness, 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 velocity-pressure 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.

2.   Governing equations

2.1.   Problem statement

Let $\Omega\subset\mathbb{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 $\mathit{\boldsymbol{u}}$, the rescaled fluid vorticity $\mathit{\boldsymbol{\omega}}$, 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. $(\mathit{\boldsymbol{x}}, t)\in\Omega_T: = \Omega\times(0, T]$, we consider

where $\mu$ is the fluid viscosity (in the considered regime, it is assumed independent of the bacteria and chemical concentrations), $\mathbb{K}(\mathit{\boldsymbol{x}})$ is the permeability tensor rescaled with viscosity, $s\mathit{\boldsymbol{g}}$ represents the force exerted by the bacteria on the fluid motion (where $\mathit{\boldsymbol{g}}$ is the gravitational acceleration, or any other particular force to which the species comply), and $\mathit{\boldsymbol{f}}(\mathit{\boldsymbol{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_\partial$ together with a compatible vorticity trace $\mathit{\boldsymbol{\omega}}_\partial$ are imposed along the domain boundary $\partial\Omega$. The presentation will be restricted however to homogeneous boundary conditions. The analysis readily extends to the non-homogeneous case by classical lifting arguments.

2.2.   Main assumptions

We suppose that the permeability tensor $\mathbb{K}\in [C(\bar{\Omega})]^{3\times 3}$ is symmetric and uniformly positive definite. So is its inverse, i.e., there exists $C>0$ such that

In addition, the diffusivities are assumed always positive, coercive, and continuous: for $i\in\{s, c\}$,

Regarding the reaction terms $G_c, G_s$, we assume they are continuous functions and there exist constants $C_c, C_s, C>0$ such that

for some parameters $\nu_1, \nu_2>0$. Initial data are assumed nonnegative and regular enough

Observe that for constant coefficients $D_c, D_s, G_c, G_s$, problem (2.1) might also serve as a model for the coupling of Newtonian flows with mass and heat transport [32].

2.3.   Weak solutions

According to (2.2), let us introduce the functional spaces

where $\mathbf{H}({\rm{div}};\Omega)$ is the space of square integrable vectorial functions whose divergence is a scalar square integrable function in $\Omega$. We then proceed, for a given $t>0$, to test the first two equations in (2.1) with functions in ${\rm{H}}^1_0(\Omega)$, the third equation in (2.1) against functions in $\mathbf{H}_0({\rm{div}};\Omega)$, the fourth equation in (2.1) against functions in $\mathbf{H}_0({\boldsymbol{\rm{curl}}};\Omega)$ and to integrate by parts. The last equation is tested with functions in ${\rm{L}}_0^2(\Omega)$ (denoting the scalar space of square integrable functions having zero mean value in $\Omega$) and no integration by parts is applied. Consequently, and using (2.2) we can give the following definition.

Definition 2.1. In view of the properties outlined in Section 2.2, we shall say that the function $(c, s, \mathit{\boldsymbol{u}}, \mathit{\boldsymbol{\omega}}, p)$ is a weak solution to (2.1) if

and

for all $m^i\in {\rm{L}}^2(0, T;{\rm{H}}^1(\Omega))$, $\mathit{\boldsymbol{v}}\in {\boldsymbol{\rm{L}}}^2(0, T;{\boldsymbol{\rm{H}}}_0({\rm{div}};\Omega))$, $\mathit{\boldsymbol{z}}\in {\boldsymbol{\rm{L}}}^2(0, T;{\boldsymbol{\rm{H}}}_0({\boldsymbol{\rm{curl}}};\Omega))$, $q\in {{\rm{L}}}^2(0, T;{{\rm{L}}}_0^2(\Omega))$.

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.

Theorem 2.2. Assume that conditions (2.3), (2.4) and (2.5) hold. If $c_{0}, {s_{0}}\in {\rm{L}}^\infty(\Omega)$ with $ c_0 \ge 0$ and $s_0\ge 0$ a.e. in $\Omega$, then there exists a weak solution of (2.1)-(2.2) in the sense of Definition 2.1. Moreover, the solution $(c, s)$ is nonnegative and bounded in ${\rm{L}}^\infty(\Omega_T, \mathbb{R}^2)$.

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

where $ \hat{s}$ 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.

3.1.   Preliminaries

Let us recall the following abstract result (see e.g. [17,Theorem 1.3]):

Theorem 3.1. Let $(\mathcal X, \langle\cdot, \cdot\rangle_{\mathcal X})$ be a Hilbert space. Let $\mathcal A: \mathcal X\times \mathcal X \to \mathbb{R}$ be a bounded symmetric bilinear form, and let $\mathcal G: \mathcal X\to\mathbb{R}$ be a bounded functional. Assume that there exists $\bar\beta > 0$ such that

Then, there exists a unique $x \in \mathcal X$, such that

Moreover, there exists $C>0$, independent of the solution, such that

Let us also consider the kernel of the bilinear form $\displaystyle {\int_{\Omega}}\;q \;{\rm{div}}\;\mathit{\boldsymbol{u}}\;{\rm{d}}\mathit{\boldsymbol{x}}$, that is

Moreover, we endow the space ${\boldsymbol{\rm{H}}}({\boldsymbol{\rm{curl}}};\Omega)$ with the following $\mu$-dependent norm:

and recall the following inf-sup condition (cf. [17]): there exists $\beta_2>0$ only depending on $\Omega$, such that

3.2.   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 ${\rm{L}}^2(\Omega_T)$:

for $\phi\in\{c, s\}$, where ${k_m} \ge \sup \{ {\left\| {{c_0}} \right\|_{{{\rm{L}}^\infty }(\Omega )}},{\left\| {{s_0}} \right\|_{{{\rm{L}}^\infty }(\Omega )}}\} $, and $\lambda, \beta$ are defined in (3.6) and Lemma 3.4, respectively. Our idea is quite simple, for a given $\hat{s}\in \mathcal{K}_s$ we find first the velocity $\mathit{\boldsymbol{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 $\hat{s} \in \mathcal{K}_s$ and for any $t>0$.

Lemma 3.2. Assume that $\hat{s}\in \mathcal{K}_s$. Then, the variational problem

admits a unique solution $(\mathit{\boldsymbol{u}}, \mathit{\boldsymbol{\omega}}, p)\in {\boldsymbol{\rm{H}}}({\rm{div}};\Omega)\times{\boldsymbol{\rm{H}}}({\boldsymbol{\rm{curl}}};\Omega)\times{\rm{L}}_{0}^2(\Omega)$. Moreover, there exists $C>0$ independent of $\mu$ such that

Proof. First, we observe that, owing to the inf-sup condition (3.1), problem (3.3) is equivalent to: Find $(\mathit{\boldsymbol{u}}, \mathit{\boldsymbol{\omega}})\in{\rm{X}}\times{\boldsymbol{\rm{H}}}({\boldsymbol{\rm{curl}}};\Omega)$ such that

Then, the desired result follows from Theorem 3.1, repeating the argument employed in the proofs of [4,Theorem 2.2 and Corollary 2.1].

On the other hand, given a fixed velocity $\mathit{\boldsymbol{u}} \in {\boldsymbol{\rm{L}}}^2(0, T;{\boldsymbol{\rm{H}}}({\rm{div}};\Omega))$, the following result establishes the solvability of the reaction-diffusion system:

Lemma 3.3. For any $\mathit{\boldsymbol{u}} \in {\boldsymbol{\rm{L}}}^2(0, T;{\boldsymbol{\rm{H}}}({\rm{div}};\Omega))$, the system

is uniquely solvable and there exists $C>0$, depending on $\|c_0\|_{0, \Omega}$ and $\|s_0\|_{0, \Omega}$, such that

Proof. It suffices to combine assumptions (2.4), (2.6) with the general result for quasilinear parabolic problems given in [22,Section 5].

As the next step, we consider a constant $\lambda>0$ and define the auxiliary variables $(\phi_c, \phi_s)$ by setting

Then $(c, s)$ satisfies the strong form of (3.5) with diffusion and reaction terms replaced, respectively, by

We now introduce a map $\mathcal{R}:\mathcal{K}_s\to \mathcal{K}_s$ such that $\mathcal{R}(\hat{s}) = s$, where $s$ solves (3.5). The goal is to prove that such map has a fixed-point. First, let us show that $\mathcal{R}$ is a continuous mapping. Let $(\hat{s}_n)_n $ be a sequence in $\mathcal{K}_s$ and $\hat{s} \in \mathcal{K}_s$ be such that $\hat{s}_n\to \hat{s}$ in ${\rm{L}}^2(\Omega_T)$ as $n \to \infty$. Let us then define $s_n = \mathcal{R}(\hat{s}_n)$, i.e., $s_{n}$ is the solution of (3.5) associated with $\hat{s}_n$ and the solution $\mathit{\boldsymbol{u}}$ of (3.3). The objective is to show that $s_n$ converges to $\mathcal{R}(\hat{s})$ in ${\rm{L}}^2(\Omega_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 $\hat{s}\in \mathcal{K}_s$. Then

$(i)$ There exists a constant $\gamma \ge 0$ such that $0\le c_n(\mathit{\boldsymbol{x}}, t), s_n(\mathit{\boldsymbol{x}}, t)\le e^{\gamma t} k_m$, for a.e. $(\mathit{\boldsymbol{x}}, t)\in \Omega_T$, where $k_m$ is defined as in (3.2).

$(ii)$ The sequence $(c_n, s_n)_n$ is bounded in ${\rm{L}}^2(0, T;{\rm{H}}^1(\Omega, \mathbb{R}^2))\cap {\rm{L}}^\infty(0, T;{\rm{L}}^2(\Omega, \mathbb{R}^2))$.

$(iii)$ The sequence $({\partial _t} c_n, {\partial _t} s_n)_n$ is bounded in ${\rm{L}}^2(0, T;({\rm{H}}^1(\Omega, \mathbb{R}^2))^{'})$.

$(iv)$ The sequence $(c_n, s_n)_n$ is relatively compact in ${\rm{L}}^2(\Omega_T)$.

Proof. $(i)$ We replace the strong form of (3.5) by

Multiplying (3.7) by $ - c_n^ - = \frac{{{c_n} - \left| {{c_n}} \right|}}{2}$ and integrating over $\Omega$, we obtain

Now, we use that ${\rm{div}}\mathit{\boldsymbol{u}} = 0$ and (2.5) to deduce

and

According to the positivity of the second and the fourth terms in the left-hand side of (3.8), we get

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\in \mathbb{R}$ such that ${k_c} \ge {\left\| {{c_0}} \right\|_{{{\rm{L}}^\infty }(\Omega )}}$ and ${k_s} \ge {\left\| {{s_0}} \right\|_{{{\rm{L}}^\infty }(\Omega )}}$. Let us consider $t \in (0, T)$, and proceed to multiply (3.7) by

and to integrate over $\Omega$, which yields that there exists some constant $C_6>0$ such that

where $\xi_s: = (s_n-e^{(\beta-\lambda) 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 $\beta \ge \lambda\ge C_6+C_7$, an application of (3.11) yields

and exploiting that $c_0, s_0\le k$ in $\Omega$, we conclude that $c_n(\cdot, t)\le e^{(\beta-\lambda) t} k$ and $s_n(\cdot, t)\le e^{(\beta-\lambda) t} k$ in $\Omega\times (0, T)$.

$(ii)$ We multiply equation (3.7) by $c_n$ and integrate over $\Omega$ to obtain

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 $\varphi \in {\rm{L}}^2(0, T;{\rm{H}}^1(\Omega))$ 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 (ⅲ).

In summary, Lemmas 3.2, and 3.4 imply that there exists

such that, up to extracting subsequences if necessary,

and consequently, the continuity of $\mathcal{R}$ on $\mathcal{K}_s$ holds. Moreover, Lemma 3.4 indicates the boundedness of $\mathcal{R}(\mathcal{K}_s)$ within the set

Appealing to the theory of compact sets [33], the compactness of the map $\mathcal{S}\hookrightarrow{\rm{L}}^{2}(\Omega_T)$ implies that of $\mathcal{R}$. Therefore, and thanks to Schauder's fixed-point theorem, the operator $\mathcal{R}$ has a fixed point $s$. That is, there exists a solution to

for all $m^i\in {\rm{L}}^2(0, T;{\rm{H}}^1(\Omega))$, $\mathit{\boldsymbol{v}}\in {\boldsymbol{\rm{L}}}^2(0, T;{\boldsymbol{\rm{H}}}_0({\rm{div}};\Omega))$, $\mathit{\boldsymbol{z}}\in {\boldsymbol{\rm{L}}}^2(0, T;{\boldsymbol{\rm{H}}}_0({\boldsymbol{\rm{curl}}};\Omega))$, $q\in{\rm{L}}^2(0, T;{\rm{L}}_0^2(\Omega))$.

4.   Uniqueness of weak solutions

The following result completes the analysis of unique solvability for the reaction-diffusion-Brinkman system.

Theorem 4.1. Assume (2.3)-(2.5) hold, and let $(c_1, s_1, \mathit{\boldsymbol{u}}_1, \mathit{\boldsymbol{\omega}}_1, p_1), (c_2, s_2, \mathit{\boldsymbol{u}}_2, \mathit{\boldsymbol{\omega}}_2, p_2)$ be two weak solutions to (2.1)-(2.2), associated with the data $c_0 = c_{1, 0}$, $s_{0} = s_{1, 0}$, $\mathit{\boldsymbol{g}} = \mathit{\boldsymbol{g}}_1$, $\mathit{\boldsymbol{f}} = \mathit{\boldsymbol{f}}_1$ and $c_{0} = c_{2, 0}$, $s_{0} = s_{2, 0}$, $\mathit{\boldsymbol{g}} = \mathit{\boldsymbol{g}}_2$, $\mathit{\boldsymbol{f}} = \mathit{\boldsymbol{f}}_2$, respectively. Then, for any $t\in (0, T]$ there exists $C>0$ such that

where ${\mathcal{C}} = c_1-c_2$, ${\mathcal{S}} = s_1-s_2$, ${\mathcal{U}} = \mathit{\boldsymbol{u}}_1-\mathit{\boldsymbol{u}}_2$, ${\mathcal{W}} = \mathit{\boldsymbol{\omega}}_1-\mathit{\boldsymbol{\omega}}_2$, and ${\mathcal{P}} = p_1-p_2$. In particular, there exists at most one weak solution to the reaction-diffusion-Brinkman system (2.1)-(2.2).

Proof. First, we observe that the $({\mathcal{U}}, {\mathcal{W}}, {\mathcal{P}})$ satisfies

for $t\in (0, T)$ and for all $\mathit{\boldsymbol{v}}\in {\boldsymbol{\rm{L}}}^2(0, T;{\boldsymbol{\rm{H}}}_0({\rm{div}};\Omega))$, $\mathit{\boldsymbol{z}}\in {\boldsymbol{\rm{L}}}^2(0, T;{\boldsymbol{\rm{H}}}_0({\boldsymbol{\rm{curl}}};\Omega))$, $q\in{\rm{L}}^2(0, T;{\rm{L}}_0^2(\Omega))$.

After substituting $\mathit{\boldsymbol{v}} = {\mathcal{U}}$, $z = {\mathcal{W}}$, $q = {\mathcal{P}}$ in (4.2) and adding the resulting equations, we can invoke the continuous dependence on the data (3.4) to establish the a priori bound

for some constant $C_{15}>0$. Next, note that $(\mathcal{C}, \mathcal{S})$ satisfies

for $t\in (0, T)$ and for all $m^i\in {\rm{L}}^2(0, T;{\rm{H}}^1(\Omega))$ and ${\partial _t} m^i\in {\rm{L}}^2(0, T;({\rm{H}}^1(\Omega))')$ with $m^c(\cdot, T) = 0$, $i = c, s$. Now, for $t_0 \in (0, T)$ we take

where the function $\mathcal{D}$ is the so-called Kirchhoff transform (see for e.g. [25]) $\displaystyle \mathcal{D}_i(i): = \int_0^i D_i(r)\, dr$, for $i = c, s$. Using these relations in (4.4), and recalling that the function $\mathcal{D}$ is strictly increasing, we get

for some constants $C_\mathcal{D}, C_{16}>0$. Moreover,

for some constants $C_{17}, C_{18}>0$ depending on $\|c_1\|_{{\rm{L}}^\infty(\Omega_T)}$ and $\|\mathit{\boldsymbol{u}}_2\|_{{\rm{L}}^\infty(\Omega_T, \mathbb{R}^3)}$ (the latter being bounded thanks to classical maximal regularity results for Stokes equations cf. [36]). Next, we use (2.5) and the ${\rm{L}}^\infty$-bounds of $(c_1, s_1)$ and $(c_2, s_2)$, leading to

for some constant $C_{19}>0$. We proceed to collect the results (4.3) and (4.5)-(4.7), to infer that

for some constants $C_{20}, C_{21}>0$. Similarly we get for the equation of $\mathcal{S}$

for some constants $C_{22}, C_{23}, C_{24}>0$. Thus (4.8) and (4.9) imply

for some constants $C_{25}, C_{26}, C_{27}>0$. Then, an application of Gronwall's inequality to (4.10) combined with (4.3) proves (4.1).

5.   Numerical approximation

In this section, we construct a primal-mixed fully-discrete 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 $\mathcal{T}_h$ be a regular family of triangulations of $\bar\Omega$ by tetrahedra $K$ of maximum diameter $h$. Given an integer $k\ge0$ and $S\subset \mathbb{R}^3$, by $\mathcal{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\in \mathcal{T}_h(\Omega)$, the local lowest order Raviart-Thomas and edge space of Nédélec type, are defined as ${\boldsymbol{\rm{RT}}}_0(K) = \mathcal{P}_{0}(K)^3\oplus\mathcal{P}_{0}(K)\mathit{\boldsymbol{x}}, $ and ${\boldsymbol{\rm{ND}}}_1(K) = \mathcal{P}_{0}(K)^3\oplus\mathcal{P}_{0}(K)^3\times \mathit{\boldsymbol{x}}$, respectively.

Then, a Galerkin semi-discretisation associated to the formulation introduced in Definition 2.1 reads: For $t\in(0, T]$ find $(c_h(t), s_h(t), \mathit{\boldsymbol{u}}_h(t), \mathit{\boldsymbol{\omega}}_h(t), p_h(t))\in {\rm{V}}_h\times{\rm{V}}_h\times{\boldsymbol{\rm{H}}}_h\times{\boldsymbol{\rm{Z}}}_h\times{\rm{Q}}_h$ such that

for all $m_h^c, m_h^s\in {\rm{V}}_h$, $\mathit{\boldsymbol{v}}_h\in {\boldsymbol{\rm{H}}}_h$, $\mathit{\boldsymbol{z}}_h\in {\boldsymbol{\rm{Z}}}_h$, and $q_h\in{\rm{Q}}_h$.

5.2.   Euler time discretisation

Let $c_h^{0} = \mathbb{P}_{{\rm{V}}_h}(c_0)$, $s_h^{0} = \mathbb{P}_{{\rm{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 $\Delta t = T/N$: For $n\in \{1, \ldots, N\}$, find $(c^n_h, s^n_h, \mathit{\boldsymbol{u}}^n_h, \mathit{\boldsymbol{\omega}}^n_h, p^n_h)\in {\rm{V}}_h\times{\rm{V}}_h\times{\boldsymbol{\rm{H}}}_h\times{\boldsymbol{\rm{Z}}}_h\times{\rm{Q}}_h$ such that

for all $m_h^c, m_h^s\in {\rm{V}}_h$, $\mathit{\boldsymbol{v}}_h\in {\boldsymbol{\rm{H}}}_h$, $\mathit{\boldsymbol{z}}_h\in {\boldsymbol{\rm{Z}}}_h$, $q_h\in{\rm{Q}}_h$, where the forcing term in the momentum equation is discretised explicitly. We also stress that, owing to the choice (5.1), the discrete velocities $\mathit{\boldsymbol{u}}_h$ generated using either (5.2) or (5.3) are exactly divergence-free, that is, $ {\rm{div}}\;\mathit{\boldsymbol{u}}_h(t) = 0$ and ${\rm{div}}\;\mathit{\boldsymbol{u}}_h^n = 0$ in $\Omega$ (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 ${\rm{L}}^2(\Omega, \mathbb{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

to treat the finite difference discretisation of the time derivatives. Proceeding in this way, we can obtain estimates (uniform in both $h$ and $\Delta t$) for the discrete Brinkman solutions

and for the advection-reaction-diffusion system

for some constant $C>0$.

The next goal is to establish the relative compactness in ${\rm{L}}^2(\Omega_{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.

Lemma 5.1. There exists a constant $C>0$ depending on $\Omega$, $T$, $c_0$ and $s_0$ such that

for all $\mathit{\boldsymbol{r}}\in\mathbb{R}^3$ and for all $\tau\in (0, T)$, where $\Omega_{\mathit{\boldsymbol{r}}} = \{\mathit{\boldsymbol{x}} \in \Omega, \, [\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{x}}+\mathit{\boldsymbol{r}}]\subset \Omega\}$.

Proof. Let us introduce the space translates

From the ${\rm{L}}^2(0, T;{\rm{H}}^1(\Omega))-$estimate for $c_{h}$ and $s_h$, the bound

easily follows. It is then clear that the right-hand side in (5.8) vanishes as $|\mathit{\boldsymbol{r}}|\to 0$, uniformly in $h$, which yields (5.6).

Next we introduce the time translates

and notice that for all $t\in(0, T-\tau]$, these functions assume values in $V_h$ (cf. (5.1)). Therefore they can be used as test functions in the fully-discrete scheme (5.3). Moreover, the previously proved uniform bounds for $(c_{h}, s_h)$ and $(\nabla c_{h}, \nabla s_{h})$ in ${\rm{L}}^2(\Omega_{T}, \mathbb{R}^2)$ imply analogous bounds for the translates $T^\tau c_{h}, T^\tau s_{h}$ and $\nabla T^\tau c_{h}, \nabla T^\tau s_{h}$ in ${\rm{L}}^2((0, T-\tau)\times\Omega)$.

Let us now define $(\overline{c}_h, \overline{s}_h)$ as the piecewise affine in $t$ function in $W^{1, \infty}((0; T];V_h)$ interpolating the states $(c^n_h, s^n_h)_{n = 0, \dots, N} \subset V_h$ at the points $(n\Delta t)_{n = 0, \dots, N}$ (recall the interpolation scheme $\displaystyle \overline{n}_h = n_h^{n-1}+\frac{t-t_{n-1}}{\Delta t}(n_h^{n}-n_h^{n-1})$ for $n = c, s$). Then we have

for all $m_h^c, m_h^s\in {\rm{V}}_h$. We integrate (5.9) with respect to time $\sigma\in [t, t+\tau]$ (with $0<\tau <T$) and consider $m_h^c = T^\tau c_h$ and $m_h^s = T^\tau 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| \leq C\, \tau$, for some constant $C>0$. Herein we have used the Fubini theorem (recall that $\int_t^{t+\tau}\, {\rm{d}}\sigma = \tau = \int_{\sigma-\tau}^s\, dt$), the divergence-free condition for the discrete velocity, the Hölder inequality and the ${\rm{L}}^2-$bounds for $(c_h, s_h)$, $(\nabla c_{h}, \nabla s_{h})$ and $(\nabla T^h c_{h}, \nabla 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| \leq C\, \tau, $ for some constant $C>0$. Collecting these inequalities we readily get

Furthermore, the definition of ($\bar c_h, \bar 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 Ⅳ.25]), we can assert that there exists a subsequence of $(c_h, s_h, \mathit{\boldsymbol{u}}_h, \mathit{\boldsymbol{\omega}}_h, p_h)$, not re-labeled, such that, as $h\to 0$,

These convergences allow us to identify the limit $(c, s, \mathit{\boldsymbol{u}}, \mathit{\boldsymbol{\omega}}, p)$ as the weak solution of (2.1), and the convergence result is summarised as follows.

Theorem 5.2. Assume that conditions (2.3), (2.4) and (2.5) hold. If $c_{0}, s_{0}\in{\rm{L}}^2(\Omega)$, then the finite element solution $(c^n_h, s^n_h, \mathit{\boldsymbol{u}}^n_h, \mathit{\boldsymbol{\omega}}^n_h, p^n_h)$, generated by (5.3), converges along a subsequence to $(c, s, \mathit{\boldsymbol{u}}, \mathit{\boldsymbol{\omega}}, p)$ as $h, \Delta t \to 0$, where $(c, s, \mathit{\boldsymbol{u}}, \mathit{\boldsymbol{\omega}}, p)$ is a weak solution of (2.1)-(2.2) in the sense of Definition 2.1.

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 $\Delta 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

with $t_n = n\Delta t$, for $n = 1, \ldots, N$ and $C_{1}, C_{2}>0$ are constants independent of $h$ and $\Delta t$. Here $(c^n_h, s^n_h, \mathit{\boldsymbol{u}}^n_h, \mathit{\boldsymbol{\omega}}^n_h, p^n_h)$ denotes the sequence generated by (5.3) for all $n = 1, \ldots, 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 reaction-diffusion 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 strategies and thorough comparisons between dedicated solvers have been reported in [24].

Example 1. Spatio-temporal accuracy. Let us consider the spatio-temporal domain $\Omega_T = (0, 1)^2\times(0, 0.1]$ and define the following exact solutions to (2.1) in the case of constant diffusivities $D_c, D_s$:

The model parameters are set as $\mu = 1$, $D_c = 0.01$, $D_s = 0.5$, $\mathit{\boldsymbol{g}} = (0, -1)^T$, and $\mathbb{K} = {\boldsymbol{\rm{I}}}$; the reaction kinetics are specified by

whereas $\mathit{\boldsymbol{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.

Example 1. Average number of Newton iterations to reach the residual tolerance tol$ = 1e-6$, according to the number of degrees of freedom in each refinement level.

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 $\Delta t = $1.6e-4 and computing errors between the exact solutions from (6.1) and approximations obtained with our mixed-primal method on a sequence of successively refined structured triangulations of $\Omega$, at the final adimensional time $T = $5e-3. That is, for each field $\eta$, we compute $ e_h(\eta): = \|\eta^N-\eta^N_h\|$, where $\eta^N = \eta(\cdot, T)$, and $\|\cdot\|$ is either the ${\rm{H}}^1$, ${\rm{H}}({\rm{div}})$, or ${\rm{L}}^2-$norm. On the other hand, the error history associated to the time discretisation is studied by fixing a small meshsize $h = 2.43$e-5 and computing accumulative errors up to $T = 0.1$, and decreasing $\Delta t$. That is, we measure errors in the $\ell^\infty(0, t;\cdot)-$norm: $E_{\Delta t}(\eta): = \sum_{n = 0}^{N}\|\eta^n-\eta_h^n\|_{(\cdot)}$. 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) = -\frac{\alpha}{2}(1+ \frac{s -s^*}{\sqrt{(s-s^*)^2 + \varepsilon^2}})$, 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) = \beta cr(s)$, and the remaining model functions and adimensional parameters are set as follows: $\mathbb{K} = K_0 + K_0 \eta(x)$ (with $K_0 = 7700$ and $\eta$ a uniformly distributed random field), $g = (0, -\gamma)^T$, $s^* = 0.3$, $D_c = 0.01$, $D_s = 0.25$, $\varepsilon = h$. The computational domain is a disk of radius $\frac{1}{2}$ centred at $(\frac{1}{2}, \frac{1}{2})$, discretised into an unstructured mesh of 13972 points and 27942 triangular elements, and a timestep of $\Delta 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\sqrt{(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 $\alpha, \beta, \gamma$. 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.

Example 3. FitzHugh-Nagumo dynamics in a porous cavity. This test emphasises the effects of high contrast permeabilities on the propagation of travelling waves dictated by the well-known FitzHugh-Nagumo reactions $ G_c(c, s) = k(s+c(c-a)(c-1))$, and $G_s(c, s) = d_1c-s$. In this prototype model of excitable systems, $c$ and $s$ represent a membrane voltage, and a recovery variable. We consider a rectangular domain $\Omega = (0, 1.8)\times(0, 1)$ filled with 60 large particles (of radii 0.0015 and located randomly). The forcing terms are $\mathit{\boldsymbol{g}} = \sqrt{2} d_2(1, 1)^T$ and $\mathit{\boldsymbol{f}} = \mathit{\boldsymbol{0}}$. Moreover, the conductivity of the medium is assumed anisotropic and affected by the heterogeneity of the permeability field $(D_c)_{1, 1} = 40\mathbb{K}$, $(D_c)_{2, 2} = 5\mathbb{K}$, $D_s = $1e-3, the fluid viscosity is $\mu = 0.01$, and the remaining parameters are chosen as $a = 0.25$, $d_1 = 0.16875$, $d_2 = 250\sqrt{3}$, $k = -100$. The mesh has 28712 elements sharing 14581 nodes, we set a timestep of $\Delta 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 $\Omega = (0, 1.8)\times(0, 1)\times(0, 0.6)$ filled with 60 large particles (of radii 1.5e-3 and randomly distributed on $\Omega$) with a permeability 100 times higher than in the rest of the porous matrix. The forcing term is specified by $\mathit{\boldsymbol{g}} = d_2(1, 1, 1)^T$ and $\mathit{\boldsymbol{f}} = \mathit{\boldsymbol{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 $\Delta 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 species-dependent coefficients:

where $\eta$ 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 $ \mathit{\boldsymbol{f}}(s) = \gamma_0 |s| (\mathit{\boldsymbol{f}}_0\otimes\mathit{\boldsymbol{f}}_0) \nabla s$, which might be regarded as the flow-counterpart to the active force proposed in [28]. Here $\mathit{\boldsymbol{f}}_0 = (1, 0)^T$ is a given preferential direction of plasma displacement, and the remaining parameters are $\gamma_0 = -0.12$, $\nu_1 = 1.58$, $\nu_2 = 16$, $\nu_3 = 91$, $\nu_4 = 2$, $\nu_5 = 0.2$, $k_2 = 1$, $k_3 = 4$, $k_4 = 0.75$, $\mathbb{K} = 1\text{e-4}(1 + 10 \eta(\mathit{\boldsymbol{x}}))$. Snapshots of the numerical solutions at an early and advanced time steps are collected in Figure 5. Starting from the open channels, the cytosolic calcium starts to propagate towards other channels, also following a higher diffusion in the preferential direction $\mathit{\boldsymbol{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.

Acknowledgments

We thank the inputs from two anonymous referees, which have resulted in a number of improvements with respect to the initial version of the manuscript. Also, special thanks are due to Davin Lunz, Tori Pereira, Lindon Roberts, Caoimhe Rooney, Ian Roper, Valentin Sulzer, and Jessica Williams (cohort two of the Oxford Centre for Doctoral Training on Industrially Focused Mathematical Modelling) for the design and preliminary implementation of Example 2 (bacterial bioconvection on a disk). This work was partially supported by CONICYT-Chile through FONDECYT projects 11160706 and 1140791; by DIUBB through project 165608--3/R; and by the EPSRC through the Research Grant EP/R00207X/1.