A virtual element method for the miscible displacement of incompressible fluids in porous media

TL;DR

This paper develops a virtual element method for simulating the displacement of incompressible fluids in porous media, providing a new numerical approach with theoretical analysis and practical testing.

Contribution

It introduces the first application of virtual element methods to complex fluid flow problems involving miscible displacement, combining discretization, time stepping, and analysis.

Findings

The method is effective on regular test cases.

The method performs well on realistic scenarios.

Theoretical analysis confirms stability and convergence.

Abstract

In the present contribution, we construct a virtual element (VE) discretization for the problem of miscible displacement of one incompressible fluid by another, described by a time-dependent coupled system of nonlinear partial differential equations. Our work represents a first study to investigate the premises of virtual element methods (VEM) for complex fluid flow problems. We combine the VEM discretization with a time stepping scheme and develop a complete theoretical analysis of the method under the assumption of a regular solution. The scheme is then tested both on a regular and on a more realistic test case.

Equations

$$\{\begin{array}{rl}\varphi \frac{\partial c}{\partial t}+\mathbf{u}\cdot \nabla c-\mathrm{div}(D(\mathbf{u})\nabla c)& =q^{+}(\hat{c}-c)\\ \mathrm{div}\mathbf{u}& =G\\ \mathbf{u}& =-a(c)(\nabla p-\mathbf{\gamma }(c)),\end{array}$$

$$G:=q^{+}-q^{-}.$$

$$D(\mathbf{u}):=\varphi \left[d_{m}I+|\mathbf{u}|(d_{\ell }E(\mathbf{u})+d_t{E}^{\perp }(\mathbf{u}))\right],$$

$$E(\mathbf{u}):={\left(\frac{{\mathbf{u}}_i{\mathbf{u}}_j}{|\mathbf{u}{|}^2}\right)}_{i,j=1,2}=\frac{\mathbf{u}{\mathbf{u}}^T}{|\mathbf{u}{|}^2},E^{\perp }(\mathbf{u}):=I-E(\mathbf{u}),$$

$$a(c):=\frac{k}{\mu (c)},$$

$$\mu (c)=\mu (0){\left(1+\left(M^{\frac{1}{4}}-1\right)c\right)}^{-4},\text{ in }[0,1],$$

$$\{\begin{array}{rl}\mathbf{u}\cdot \mathbf{n}& =0\text{ on }\partial \Omega \times J\\ D(\mathbf{u})\nabla c\cdot \mathbf{n}& =0\text{ on }\partial \Omega \times J,\end{array}$$

$$c(\mathbf{x},0)=c_0(\mathbf{x})\text{ in }\Omega ,$$

$${\int }_{\Omega }{q}^{+}(\mathbf{x},t)\text{d}x={\int }_{\Omega }{q}^{-}(\mathbf{x},t)\text{d}x,$$

$$a_{\ast }\le a(z,\mathbf{x})\le a^{\ast },{\varphi }_{\ast }\le \varphi (\mathbf{x})\le {\varphi }^{\ast },$$

$$A(z)(\mathbf{x}):=a^{-1}(z,\mathbf{x}).$$

$$0<d_m\le d_t\le d_{\ell }.$$

$$H(\mathrm{div};\Omega ):=\{\mathbf{v}\in [L^2(\Omega ){]}^2:\mathrm{div}\mathbf{v}\in L^2(\Omega )\}.$$

$$\begin{array}{rl}\mathbf{V}& :=\{\mathbf{v}\in H(\mathrm{div};\Omega ):\mathbf{v}\cdot \mathbf{n}=0\text{ on }\partial \Omega \}\\ Q& :=L_0^2(\Omega ):=\{\phi \in L^2(\Omega ):(\phi ,1{)}_{0,\Omega }=0\}\\ Z& :=H^1(\Omega ),\end{array}$$

$$\begin{array}{r}\Vert \mathbf{u}{\Vert }_{\mathbf{V}}^2:=\Vert \mathbf{u}{\Vert }_{0,\Omega }^2+\Vert \mathrm{div}\mathbf{u}{\Vert }_{0,\Omega }^2,\Vert q{\Vert }_Q^2:=\Vert q{\Vert }_{0,\Omega }^2,\Vert z{\Vert }_Z^2:=\Vert z{\Vert }_{1,\Omega }^2:=\Vert \nabla z{\Vert }_{0,\Omega }^2+\Vert z{\Vert }_{0,\Omega }^2.\end{array}$$

$$\mathbf{u}(t)(x):=\mathbf{u}(x,t),p(t)(x):=p(x,t),c(t)(x):=c(x,t).$$

$$\Vert \mathbf{v}{\Vert }_{L^2(a,b;\mathbf{V})}:={\left({\int }_a^b\Vert \mathbf{v}(t){\Vert }_{\mathbf{V}}^2\text{d}x\right)}^{\frac{1}{2}},\Vert \mathbf{v}{\Vert }_{L^{\infty }(a,b;\mathbf{V})}:=\underset{t\in [a,b]}{\text{ess sup}}\Vert \mathbf{v}(t){\Vert }_{\mathbf{V}};$$

$$\{\begin{array}{rl}\mathcal{M}\left(\frac{\partial c(t)}{\partial t},z\right)+{\left(\mathbf{u}(t)\cdot \nabla c(t),z\right)}_{0,\Omega }+\mathcal{D}(\mathbf{u}(t);c(t),z)& ={\left(q^{+}(\hat{c}-c)(t),z\right)}_{0,\Omega }\\ \mathcal{A}(c(t);\mathbf{u}(t),\mathbf{v})+B(\mathbf{v},p(t))& =(\mathbf{\gamma }(c(t)),\mathbf{v}{)}_{0,\Omega }\\ B(\mathbf{u}(t),q)& =-{\left(G(t),q\right)}_{0,\Omega }\end{array}$$

$$\mathcal{M}(c,z)$$

$$\mathcal{A}(c;\mathbf{u},\mathbf{v})$$

$$\begin{array}{r}\mathcal{M}\left(\frac{\partial c}{\partial t},z\right)+\Theta (\mathbf{u},c;z)+\mathcal{D}(\mathbf{u};c,z)={\left(q^{+}\hat{c},z\right)}_{0,\Omega },\end{array}$$

$$\Theta(\boldsymbol{u},c;z):=\frac{1}{2}\bigg{[}\left(\boldsymbol{u}\cdot\nabla c,z\right)_{0,\Omega}+((q^{+}+q^{-})\,c,z)_{0,\Omega}-\left(\boldsymbol{u},c\,\nabla z\right)_{0,\Omega}\bigg{]}.$$

$$\left(\boldsymbol{u}\cdot\nabla c,z\right)_{0,\Omega}=\frac{1}{2}\big{[}\left(\boldsymbol{u}\cdot\nabla c,z\right)_{0,\Omega}-(G,c\,z)_{0,\Omega}-\left(\boldsymbol{u},c\,\nabla z\right)_{0,\Omega}\big{]},$$

$$\mathcal{M}(c,z)\le {\varphi }^{\ast }\Vert c{\Vert }_{0,\Omega }\Vert z{\Vert }_{0,\Omega },\mathcal{M}(z,z)\ge {\varphi }_{\ast }\Vert z{\Vert }_{0,\Omega }^2,$$

$$\mathcal{A}(c;\mathbf{u},\mathbf{v})\le \frac{1}{a_{\ast }}\Vert \mathbf{u}{\Vert }_{0,\Omega }\Vert \mathbf{v}{\Vert }_{0,\Omega }.$$

$$\mathcal{A}(c;\mathbf{u},\mathbf{v})\le \Vert A(c){\Vert }_{0,\Omega }\Vert \mathbf{u}{\Vert }_{\infty ,\Omega }\Vert \mathbf{v}{\Vert }_{0,\Omega }.$$

$$\mathcal{A}(c;\mathbf{v},\mathbf{v})\ge \frac{1}{a^{\ast }}\Vert \mathbf{v}{\Vert }_{0,\Omega }^2$$

$$\mathcal{K}:=\{\mathbf{v}\in \mathbf{V}:B(\mathbf{v},q)=0\forall q\in Q\},$$

$$\begin{array}{r}\mathcal{D}(\mathbf{u};c,z)\le {\varphi }^{\ast }\left[d_m+\Vert \mathbf{u}{\Vert }_{\infty ,\Omega }(d_{\ell }+d_t)\right]\Vert \nabla c{\Vert }_{0,\Omega }\Vert \nabla z{\Vert }_{0,\Omega },\end{array}$$

$$\mathcal{D}(\mathbf{u};c,z)\le \Vert D(\mathbf{u}){\Vert }_{0,\Omega }\Vert \nabla c{\Vert }_{\infty ,\Omega }\Vert \nabla z{\Vert }_{0,\Omega }\le {\eta }_{\mathcal{D}}(1+\Vert \mathbf{u}{\Vert }_{0,\Omega })\Vert \nabla c{\Vert }_{\infty ,\Omega }\Vert \nabla z{\Vert }_{0,\Omega },$$

$$\begin{array}{rl}(D(\mathbf{u})\mathbf{\mu },\mathbf{\mu }{)}_{0,\Omega }& =(\varphi d_m\mathbf{\mu },\mathbf{\mu }{)}_{0,\Omega }+(\varphi |\mathbf{u}|(d_{\ell }E(\mathbf{u})+d_t{E}^{\perp }(\mathbf{u}))\mathbf{\mu },\mathbf{\mu }{)}_{0,\Omega }\\ & \ge {\varphi }_{\ast }{d}_m\Vert \mathbf{\mu }{\Vert }_{0,\Omega }^2+(\varphi |\mathbf{u}|(d_{\ell }-d_t)E(\mathbf{u})\mathbf{\mu },\mathbf{\mu }{)}_{0,\Omega }+(\varphi |\mathbf{u}|d_t\mathbf{\mu },\mathbf{\mu }{)}_{0,\Omega }\\ & \ge {\varphi }_{\ast }\left(d_m\Vert \mathbf{\mu }{\Vert }_{0,\Omega }^2+d_t\Vert |\mathbf{u}{|}^{\frac{1}{2}}\mathbf{\mu }{\Vert }_{0,\Omega }^2\right)\end{array}$$

Full text

AMS subject classification: 65M12, 65M60, 76S05

Keywords: virtual element methods, miscible fluid flow, porous media, polygonal meshes

1 Introduction

The virtual element method (VEM) was introduced in [5, 7] (see also [14]) as a generalization of the finite element method (FEM) that allows to use general polygonal and polyhedral meshes. Since its recent birth in 2013, VEM enjoyed a rapid growth in the mathematics and engineering communities. Among the large numbers of papers in the literature, we here cite only [28, 3, 32, 22, 10, 29, 6, 20, 51, 11, 50, 45, 36] as representatives, see also the references therein.

In the realm of diffusion problems, virtual elements have been developed for linear model diffusion-convection-reaction equations in primal and mixed form [5, 9, 8, 14, 31, 58, 17]. It was soon recognized that the flexibility of VEM in terms of meshing could lead to appealing advantages in the presence of complex geometries, such as for discrete fracture network simulation [18, 16, 21] and, more in general, in the presence of fractures in porous 3D media [19, 43]. Nevertheless, although in other frameworks (such as solid mechanics) VEM have indeed proven themselves also on tough nonlinear problems, to the best knowledge of the authors, virtual elements have never been developed and tested for more complex diffusion models. Since VEM have indeed been hardly tested (with very promising outcomes) on linear diffusion problems with complex geometries, as often encountered in geophysical flows, developing a VEM also for more complex (and realistic) flow models becomes a key step towards a competitive methodology for applications.

In the present contribution, we consider the miscible displacement of one incompressible fluid by another in a reservoir, described by a time-dependent coupled system of nonlinear partial differential equations, that is a basic (but meaningful) model instrumental to applications such as oil recovery and environmental pollution [55, 52, 37, 33, 38]. One must note that, on this and similar models, there already exists a large literature with many competitive schemes, adopting for instance finite elements [40, 41, 37], discontinuous Galerkin methods [4, 54, 49], and finite volumes [33, 2, 34]. The aim of the present paper is to make a first study on the premises of VEM in this framework (by proposing a numerical scheme, giving a theoretical backbone to it, and finally testing it numerically). We believe that the VEM for this kind of problems could have a value due to its strong robustness with respect to the mesh features and, more heuristically, its potential in terms of flexibility in general terms.

From the mathematical viewpoint, the above model yields a nonlinear time-dependent coupled problem for concentration, velocity and pressure, also with potential issues of stability (at the discrete level) due to possible convection-dominated regimes. We propose a continuous ($H^{1}$ conforming) approximation for the concentration variable, thus leading to nodal virtual elements [5, 9], and an $H({\rm div})$ conforming approximation of the Darcy velocity, thus leading to face virtual elements [27, 8]. For the pressure, we adopt a standard piecewise discontinuous polynomial space. Due to the presence of the non-linear coefficients coupling the two set of equations, we make use of projection operators to approximate such terms and of stabilization factors that are suitably chosen. We combine the VEM discretization in space with a simple discretization procedure in time, that is a backward Euler approximation that is explicit in the coefficient terms. As a consequence, the system to be solved at each time step is linear and decoupled, leading to a cheap procedure. Extending the proposed scheme to different time discretization procedures would be, on the basis of the work presented here, quite trivial.

After proposing the method, we develop an error analysis under the assumption of a regular solution. Although such regularity conditions are unrealistic in most cases of interest, we believe that the derived results are still critical in order to give a theoretical backbone to the method. They serve the purpose of showing that the method indeed delivers a solution with the potential to yield the correct approximation order whenever this is feasible (given the approximability of the target solution by the discrete space). No time step size condition is needed in the analysis. Finally, we test the proposed scheme in two different ways. We firstly consider a problem with known regular solution inspired from [44], in order to validate the convergence properties of the method also in practice and to test some other practical aspect such as the possibility of having different time step sizes for the two different equations. Then, we consider a more realistic test, taken from [59], in order to have a qualitative comparison with the expected benchmark solution from the literature. In this second test, there is also the risk of overshoots and undershoots in the discrete solution due to strong convection. We here deal with this aspect by introducing in our scheme a modification borrowed from [47], that is recognized in [46] to be one of the best choices in practice. From the present first theoretical and numerical studies, we believe the VEM is promising and has the possibility to become, after further developments, a competitive scheme for complex flow problems.

The structure of the paper is the following. In Section 2, we introduce the continuous problem, in strong and weak form. In Section 3, we describe the proposed virtual element discretization, in space and time. In Section 4, we develop the theoretical convergence analysis of the scheme. Finally, in Section 5, we show the numerical results.

2 Problem description

We consider the miscible displacement of one incompressible fluid by another in a porous medium. This problem can be formulated in terms of a system of partial differential equations, where a parabolic diffusion-convection-reaction type equation is nonlinearly coupled with an elliptic system, see also [52, 37, 33, 38].

We need to introduce some notation and conventions to be adopted throughout the paper. We denote by $\mathbb{N}$ and $\mathbb{N}_{0}$ the sets of all natural numbers without and including zero, respectively. Moreover, we employ the standard notation for Sobolev spaces, norms, and seminorms. More precisely, for a given bounded Lipschitz domain $D\subset\mathbb{R}^{2}$, $k\in\mathbb{N}\cup\{\infty\}$, and $p\in\mathbb{N}$, we define by $W^{k,p}(D)$ the space of all $L^{p}$ integrable functions over $D$ whose weak derivatives up to order $k$ are again $L^{p}$ integrable. Sobolev spaces with fractional order can be defined for instance via interpolation theory [57]. For $p=2$, we write $H^{k}(D):=W^{k,2}(D)$, and we use $(\cdot,\cdot)_{k,D}$, $|\cdot|_{k,D}$, and $\lVert\cdot\rVert_{k,D}$, to denote the corresponding inner product, seminorm, and norm, respectively. The standard $L^{2}$ inner product over $D$ is written as $(\cdot,\cdot)_{0,D}$ with corresponding norm $\lVert\cdot\rVert_{0,\Omega}$. Further, $\mathbb{P}_{k}(D)$ is the space of polynomials up to order $k$, and $[\mathbb{P}_{k}(D)]^{2}$ the corresponding vector valued space. Additionally, $|\cdot|$ is the standard Euclidean norm for scalars and vectors. Finally, throughout the paper, $\eta$ denotes a generic constant, possibly varying from one occurrence to the other, but independent of the mesh size and, apart from Theorem 4.8, also independent of the variables.

2.1 Continuous Problem

Let $\Omega\subset\mathbb{R}^{2}$ be a polygonal bounded, convex domain, describing a reservoir of unit thickness. Given a time interval $J:=[0,T]$, for $T>0$, we are interested in finding $\boldsymbol{u}=\boldsymbol{u}(\boldsymbol{x},t)$, representing the Darcy velocity (volume of fluid flowing cross a unit across-section per unit time), the pressure $p=p(\boldsymbol{x},t)$ in the fluid mixture, and the concentration $c=c(\boldsymbol{x},t)$ of one of the fluids (amount of the fluid per unit volume in the fluid mixture), with $(\boldsymbol{x},t)\in\Omega_{T}:=\Omega\times J$, such that

[TABLE]

where $\phi=\phi(\boldsymbol{x})$ is the porosity of the medium, $q^{+}=q^{+}(\boldsymbol{x},t)$ and $q^{-}=q^{-}(\boldsymbol{x},t)$ are the (non negative) injection and production source terms, respectively, $\widehat{c}=\widehat{c}(\boldsymbol{x},t)$ is the concentration of the injected fluid, and

[TABLE]

Moreover, $D(\boldsymbol{u})\in\mathbb{R}^{2\times 2}$ is the diffusion tensor given by

[TABLE]

with matrices

[TABLE]

and molecular diffusion coefficient $d_{m}$, longitudinal dispersion coefficient $d_{\ell}$, and transversal dispersion coefficient $d_{t}$. Further, $\boldsymbol{\gamma}(c)$ in (1) describes the force density due to gravity (typically written as $\boldsymbol{\gamma}(c)=\gamma_{0}(c)\boldsymbol{\rho}$ with $\gamma_{0}(c)$ being the density of the fluid and $\boldsymbol{\rho}$ the gravitational acceleration vector), and $a(c)=a(c,\boldsymbol{x})$ is the scalar valued function given by

[TABLE]

where $k=k(\boldsymbol{x})$ represents the permeability of the porous rock, and $\mu(c)$ is the viscosity of the fluid mixture, which can be modeled by

[TABLE]

with mobility ratio $M:=\frac{\mu(0)}{\mu(1)}$. Note that $\mu$ can be set to $\mu(0)$ for $c<0$, and to $\mu(1)$ for $c>1$. We also highlight that, in the literature, $k$ is sometimes assumed to be a tensor. The following analysis can be straightforwardly generalized to that case.

Assuming impermeability of $\partial\Omega$, the system (1) is closed by requiring no-flow boundary conditions of the form

[TABLE]

and initial condition

[TABLE]

where $0\leq c_{0}(\boldsymbol{x})\leq 1$ is an initial concentration.

By use of the divergence theorem, the boundary conditions (4) directly imply the following compatibility condition for $q^{+}$ and $q^{-}$:

[TABLE]

for every $t\in J$.

We highlight that, in the forthcoming theoretical analysis, we will always assume sufficient regularity of the exact solution and the involved functions, such as $q^{+}$, $q^{-}$, $\widehat{c}$, et cetera, as better motivated in the corresponding section. Moreover, we will make use of the following assumptions.

First of all, we suppose that the functions $a$ and $\phi$ are positive and uniformly bounded from below and above, i.e. there exist positive constants $a_{\ast}$, $a^{\ast}$, $\phi_{\ast}$, and $\phi^{\ast}$, such that

[TABLE]

for all $\boldsymbol{x}\in\Omega$ and $z=z(t)$. For the sake of readability, we define

[TABLE]

Additionally, we will make use of the following relation of the diffusion and dispersion coefficients, which was observed in laboratory experiments:

[TABLE]

Finally, we recall that the source terms $q^{+}$ and $q^{-}$ are, as usual, assumed to be non-negative functions.

Existence of weak solutions to this model problem was shown in [42] for $\boldsymbol{\gamma}(c)=0$. An extension of this result to 3D spatial domains, including the presence of $\boldsymbol{\gamma}(c)$ and various boundary conditions was discussed in [35].

2.2 Weak formulation of the continuous problem

Here, we fix the basic notation and the functional setting.

To this purpose, given $\Omega$ as above, we first introduce the Sobolev space

[TABLE]

Then, we define the velocity space $\boldsymbol{V}$, the pressure space $Q$, and the concentration space $Z$ by

[TABLE]

respectively. These spaces are endowed, respectively, with the following norms:

[TABLE]

Note that $\operatorname{\rm div}\boldsymbol{V}=Q$.

As usual in the framework of parabolic problems, we use the notation

[TABLE]

For $0\leq a\leq b$, we further introduce

[TABLE]

analogously for $p$ and $c$.

Having this, the continuous problem reads as follows: find $c\in L^{2}(0,T;Z)\cap C^{0}([0,T];L^{2}(\Omega))$, $\boldsymbol{u}\in L^{2}(0,T;\boldsymbol{V})$, and $p\in L^{2}(0,T;Q)$, such that

[TABLE]

for all $\boldsymbol{v}\in\boldsymbol{V}$, $q\in Q$, and $z\in Z$, for almost all $t\in J$ and with initial condition $c(0)=c_{0}$, where

[TABLE]

Note that $c\in L^{2}(0,T;Z)\cap C^{0}([0,T];L^{2}(\Omega))$ implies $\frac{\partial c}{\partial t}\in L^{2}(0,T;Z^{\prime})$, see e.g. [53, Thm. 11.1.1].

For the sake of readability, we suppressed $(t)$ in (11). From now on, we will use the convention that by writing $\boldsymbol{u}$, we mean in fact $\boldsymbol{u}(t)$; similarly for the other functions depending on space and time. In general it will be clear from the context whether $\boldsymbol{u}$ represents $\boldsymbol{u}(t)$ for a fixed $t\in J$, i.e. as a function of space only, or for varying $\boldsymbol{x}$ and $t$, as a function of both space and time.

Moreover, we will use the following alternative form for the concentration equation:

[TABLE]

where

[TABLE]

This version is obtained from the original one in (10) by rewriting the convective term as

[TABLE]

where we first integrated by parts, then employed the fact that $\nabla\cdot\boldsymbol{u}=G$, together with the definition of $G$ in (2), and afterwards combined this term with $(q^{+}\,c,z)_{0,\Omega}$ from the right hand side of (10). This representation was inspired by the theory of VEM for general elliptic problems [32] and helps to ensure that properties of the continuous bilinear will be preserved after discretization.

In the rest of this section, we summarize some properties of the forms $\mathcal{M}(\cdot,\cdot)$, $\mathcal{A}(\cdot,\cdot,\cdot)$ and $\mathcal{D}(\cdot;\cdot,\cdot)$, all defined in (11), which will be needed later on.

To start with, for $\mathcal{M}(\cdot,\cdot)$, it directly holds with the Cauchy-Schwarz inequality and (6)

[TABLE]

for all $c,z\in Z$.

Concerning $\mathcal{A}(\cdot;\cdot,\cdot)$, again employing (6), for all $c\in L^{\infty}(\Omega)$ and $\boldsymbol{u},\boldsymbol{v}\in[L^{2}(\Omega)]^{2}$, we have

[TABLE]

Further, if $c\in L^{2}(\Omega)$, $\boldsymbol{u}\in[L^{\infty}(\Omega)]^{2}$ and $\boldsymbol{v}\in[L^{2}(\Omega)]^{2}$, it holds true that

[TABLE]

We also have the coercivity bound

[TABLE]

for all $c\in L^{\infty}(\Omega)$ and $\boldsymbol{v}\in[L^{2}(\Omega)]^{2}$, from which, after defining the kernel

[TABLE]

coercivity of $\mathcal{A}(c;\cdot,\cdot)$ on $\mathcal{K}$ in the norm $\lVert\cdot\rVert_{V}$ follows.

Regarding $\mathcal{D}(\cdot;\cdot,\cdot)$, the following continuity properties can be shown. Firstly, for all $\boldsymbol{u}\in[L^{\infty}(\Omega)]^{2}$ and $c,z\in H^{1}(\Omega)$, we have

[TABLE]

which follows directly from the Cauchy-Schwarz inequality, the definition of $D(\boldsymbol{u})$ in (3), and the fact that $|E(\boldsymbol{u})\boldsymbol{v}|\leq|\boldsymbol{v}|$ and $|E^{\perp}(\boldsymbol{u})\boldsymbol{v}|\leq|\boldsymbol{v}|$ for all $\boldsymbol{v}\in\mathbb{R}^{2}$. Moreover, for all $\boldsymbol{u}\in[L^{2}(\Omega)]^{2}$ and $c,z\in H^{1}(\Omega)$ with $\nabla c\in L^{\infty}(\Omega)$, we have the bound

[TABLE]

with matrix norm $\lVert D(\boldsymbol{u})\rVert_{0,\Omega}:=\left(\sum_{i,j=1}^{2}\lVert D_{i,j}(\boldsymbol{u})\rVert^{2}_{0,\Omega}\right)^{\frac{1}{2}}$, and some positive constant $\eta_{\mathcal{D}}$ depending only on $d_{m}$, $d_{\ell}$, and $d_{t}$. In addition, coercivity of $\mathcal{D}(\boldsymbol{u};\cdot,\cdot)$ for all $\boldsymbol{u}\in[L^{\infty}(\Omega)]^{2}$, with respect to $\lVert\cdot\rVert_{0,\Omega}$, follows from

[TABLE]

for all $\boldsymbol{\mu}\in[L^{2}(\Omega)]^{2}$, where we also employed (6) and (7).

3 The virtual element method

In this section, we derive a virtual element formulation for the model problem (10). To this purpose, we firstly fix the concept of polygonal decompositions of $\Omega$ in Section 3.1, and then, we introduce a set of discrete spaces, discrete bilinear forms, and projectors in Section 3.2. Having these ingredients, we state a semidiscrete formulation which is continuous in time and discrete in space in Section 3.3. The fully discrete formulation is the subject of Section 3.4.

3.1 Polygonal decompositions

Let $\mathcal{T}_{h}$ be a discretization of $\Omega$ into polygons $K$. We denote by $\mathcal{E}_{h}$ the set of all edges of $\mathcal{T}_{h}$, and, for a given element $K\in\mathcal{T}_{h}$, by $\mathcal{E}^{K}$ the set of edges belonging to $K$. Furthermore, $n_{K}$ is the number of edges of $K$, $h_{K}$ is the diameter of $K$, and $h:=\max_{K\in\mathcal{T}_{h}}h_{K}$. For a given edge $e\in\mathcal{E}_{h}$, we write $h_{e}$ for its length. Having this, we make the following assumptions on $\mathcal{T}_{h}$: there exists $\rho_{0}>0$ such that, for all $h>0$ and for all $K\in\mathcal{T}_{h}$,

• (D1)

$K$ is star-shaped with respect to a ball of radius $\rho\geq\rho_{0}h_{K}$;

• (D2)

$h_{e}\geq\rho_{0}h_{K}$ for all $e\in\mathcal{E}^{K}$.

Note that these two assumptions imply that the number of edges of each element is uniformly bounded. Additionally, we will require quasi-uniformity:

• (D3)

for all $h>0$ and for all $K\in\mathcal{T}_{h}$, it holds $h_{K}\geq\rho_{1}h$, for some positive uniform constant $\rho_{1}$.

Given $\mathcal{T}_{h}$, we define, for all $s>0$, the broken Sobolev spaces on $\mathcal{T}_{h}$ as

[TABLE]

together with the corresponding broken seminorms and norms

[TABLE]

Remark 1*.*

Both assumptions (D1) and (D2) are standard in the virtual element literature. While condition (D1) is quite critical in the following analysis, assumption (D2) could be possiby avoided by following steps similar to [15, 25], at the expense of making the proofs even more lenghty and technical. Finally, assumption (D3) (that can be found also in many FEM papers on the same subject) is only needed to prove bound (61) below.

3.2 Discrete spaces and projectors

Here, we introduce the local discrete VE spaces corresponding to $\boldsymbol{V}$, $Q$ and $Z$ in (8), a set of local projectors mapping from these VE spaces into spaces made of polynomials, and finally, the related global counterparts.

3.2.1 Local discrete spaces

Let $K\in\mathcal{T}_{h}$ and let $k\in\mathbb{N}_{0}$ be a given degree of accuracy. Then, the local velocity and pressure VE spaces are defined by

[TABLE]

These spaces are coupled with the preliminary local concentration space

[TABLE]

Moreover, it is important to observe that $[\mathbb{P}_{k}(K)]^{2}\subset\boldsymbol{V}_{h}(K)$ and $\mathbb{P}_{k+1}(K)\subseteq\widetilde{Z_{h}}(K)$. Associated sets of local degrees of freedom are given as follows:

• •

for $\boldsymbol{V}_{h}(K)$, a set of degrees of freedom $\{\textup{dof}^{\boldsymbol{V}_{h}(K)}_{j}\}_{j=1}^{{\textrm{dim}{V}_{h}(K)}}$ is defined by

[TABLE]

with $\boldsymbol{x}^{\perp}:=(\boldsymbol{x}_{2},-\boldsymbol{x}_{1})^{T}$, where we assume the coordinates to be centered at the barycenter of the element;

• •

for $Q_{h}(K)$, we consider $\{\textup{dof}^{Q_{h}(K)}_{j}\}_{j=1}^{\textrm{dim}{Q_{h}(K)}}$ with

[TABLE]

• •

for $\widetilde{Z_{h}}(K)$, we take $\{\textup{dof}^{\widetilde{Z_{h}}(K)}_{j}\}_{j=1}^{\textrm{dim}{\widetilde{Z_{h}}(K)}}$ with

[TABLE]

In all three cases, unisolvency is provided. More precisely, for $\boldsymbol{V}_{h}(K)$, this was proven in e.g. [12], for $Q_{h}(K)$ it is immediate, and for $\widetilde{Z_{h}}(K)$, see e.g. [5].

We also highlight that $\boldsymbol{V}_{h}(K)$ endowed with (20) mimics the Raviart-Thomas element, but in fact those two elements only coincide in the special case of triangles and $k=0$. An analogous result is true for $\widetilde{Z_{h}}(K)$, when compared to finite elements.

Remark 2*.*

We note that, for $k=0$, one obtains the lowest order local VE spaces. More precisely, in this case, the velocity space $\boldsymbol{V}_{h}(K)$ consists of all rotation free vector fields with constant divergence and edgewise constant normal traces, the pressure space $Q_{h}(K)$ only contains the constant functions, and the concentration space $\widetilde{Z_{h}}(K)$ is made of all harmonic functions that are linear on each edge. This motivates the choice of the present polynomial degrees for the spaces. However, in general, it is also possible to choose a degree of accuracy $k_{1}$ for $\boldsymbol{V}_{h}(K)$ and $Q_{h}(K)$, and another strictly positive one $k_{2}$ for $\widetilde{Z_{h}}(K)$; see e.g. [37] for FEM. The following analysis can be extended easily to such more general case just by keeping track of the different polynomial degrees.

Remark 3*.*

In order to really have a set of degrees of freedom in the computer code, one clearly needs to choose a basis for the polynomial test spaces appearing in (20) and (22). We here assume to take the classical choice, that is any monomial basis $\{m_{1},m_{2},..,m_{\ell}\}$ of the polynomial space satisfying $\|m_{i}\|_{L^{\infty}}\simeq 1$, $i=1,2,..,\ell$, where the $L^{\infty}$ norm has to be taken over the corresponding edge or bulk.

3.2.2 Local projections

For the construction of the method, we will need some tools to deal with VE functions due to the lack of their explicit knowledge in closed form. These tools will be provided in the form of local operators mapping VE functions onto polynomials. To this purpose, following [5, 7], we introduce the subsequent projectors.

The projector $\boldsymbol{\Pi^{0,K}_{k}}:\,[L^{2}(K)]^{2}\to[\mathbb{P}_{k}(K)]^{2}$ is defined as the $L^{2}$ projector onto vector valued polynomials of degree at most $k$ in each component: Given $\boldsymbol{f}\in[L^{2}(\Omega)]^{2}$,

[TABLE]

It can be shown, see [9], that this operator is computable for functions in $\boldsymbol{V}_{h}(K)$ only by knowing their values at the degrees of freedom (20). Moreover, one has computability also for functions of the form $\nabla z_{h}$ with $z_{h}\in\widetilde{Z_{h}}(K)$. This can be seen by using integration by parts:

[TABLE]

for all $\boldsymbol{p}_{k}\in[\mathbb{P}_{k}(K)]^{2}$, where the right hand side is computable by means of (22).

The projector $\Pi^{\nabla,K}_{k+1}:\,H^{1}(K)\to\mathbb{P}_{k+1}(K)$ is given, for every $z\in H^{1}(K)$, by

[TABLE]

where the second identity is needed to fix the constants. Computability of this mapping for functions in $\widetilde{Z_{h}}(K)$ was shown in [5, 7].

3.2.3 Discrete space for concentrations

The space introduced in (19) was a preliminary space, useful to introduce the main idea of the construction. Nevertheless, we will here make use of a more advanced space for the discrete concentration variable. Indeed, one can use the operator $\Pi^{\nabla,K}_{k+1}$ to pinpoint the local enhanced space

[TABLE]

where $\mathbb{P}_{k+1}/\mathbb{P}_{k-1}(K)$ is the space of polynomials in $\mathbb{P}_{k+1}(K)$ which are $L^{2}(K)$ orthogonal to $\mathbb{P}_{k-1}(K)$. It can be shown that the space $Z_{h}(K)$ has the same dimension and the same degrees of freedom (22) as $\widetilde{Z_{h}}(K)$, see [1, 13]. The advantage of the space $Z_{h}(K)$, when compared to $\widetilde{Z_{h}}(K)$, is that also the $L^{2}$ projector $\Pi^{0,K}_{k+1}:\,L^{2}(K)\to\mathbb{P}_{k+1}(K)$ onto polynomials of degree at most $k+1$, defined analogously to (23), is computable [7]

Finally, we state the following approximation result for the three projectors above [9, Lemma 5.1]:

Lemma 3.1**.**

Given $K\in\mathcal{T}_{h}$, let $\psi$ and $\boldsymbol{\psi}$ be sufficiently smooth scalar and vector valued functions, respectively. Then, it holds, for all $k\in\mathbb{N}_{0}$,

[TABLE]

where $\zeta>0$ only depends on the shape-regularity parameter $\rho_{0}$ in assumption (D1), and $k$.

3.2.4 Global discrete spaces and projectors

The global discrete spaces are defined via their local counterparts:

[TABLE]

with the obvious sets of global degrees of freedom.

In addition to the broken Sobolev norm (17), we introduce, for all $\boldsymbol{u}_{h}\in\boldsymbol{V}_{h}$,

[TABLE]

Moreover, we will denote by $\boldsymbol{\Pi^{0}_{k}}$, $\Pi^{\nabla}_{k+1}$ and $\Pi^{0}_{k+1}$, the global projectors which are defined elementwise as the corresponding local ones in Section 3.2.2 and 3.2.3.

The sets of global degrees of freedom $\{\textup{dof}^{\boldsymbol{V}_{h}}_{j}\}_{j=1}^{{\textrm{dim}{V}_{h}}}$, $\{\textup{dof}^{Q_{h}}_{j}\}_{j=1}^{\textrm{dim}{Q_{h}}}$, and $\{\textup{dof}^{Z_{h}}_{j}\}_{j=1}^{\textrm{dim}{Z_{h}}}$ are obtained by coupling the local counterparts given in (20), (21), and (22), respectively.

3.3 Semidiscrete formulation

Our aim in this section is to find a semidiscrete formulation for (10) which is continuous in time and discrete in space. To this purpose, we employ the same notation for the numerical approximants $\boldsymbol{u}_{h}$, $p_{h}$, and $c_{h}$, as in (9) for $\boldsymbol{u}$, $p$, and $c$, namely

[TABLE]

where the dependence on $(t)$ will be again suppressed in the sequel.

A semidiscrete variational formulation for (10) can then be written in an abstract way as follows: for almost every $t\in J$, find $\boldsymbol{u}_{h}\in\boldsymbol{V}_{h}$, $p_{h}\in Q_{h}$, and $c_{h}\in Z_{h}$, such that

[TABLE]

for all $\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}$, $q_{h}\in Q_{h}$, and $z_{h}\in Z_{h}$, and the initial condition

[TABLE]

is satisfied, where $I_{h}c_{0}$ is the VEM interpolant of $c_{0}$ in $Z_{h}$, and where the involved forms and terms in (24) are specified in the forthcoming lines.

Starting from the continuous problem (10), by simply replacing the continuous functions by their discrete counterparts, most of the resulting terms cannot be computed any more, owing to the fact that VE functions are not known explicitly in closed form. Thus, these terms need to be substituted by computable versions in the spirit of the VEM philosophy. To this purpose, the following replacements were made:

• •

The term $\mathcal{M}\left(\frac{\partial c_{h}}{\partial t},z_{h}\right)$ in the concentration equation was replaced by

[TABLE]

where the local contributions are given as

[TABLE]

with $S^{K}_{\mathcal{M}}(\cdot,\cdot)$ denoting a stabilization term with certain properties and a constant $\nu_{\mathcal{M}}^{K}(\phi)$, both described in Section 3.3.1 below.

• •

Next, the term $\Theta(\boldsymbol{u}_{h},c_{h};z_{h})$ was substituted by

[TABLE]

where

[TABLE]

• •

Moreover, the term $\mathcal{D}(\boldsymbol{u}_{h};c_{h},z_{h})$ was replaced by

[TABLE]

with local contributions

[TABLE]

where $S^{K}_{D}(\cdot,\cdot)$ is a stabilization term with certain properties and a constant $\nu_{D}^{K}(\boldsymbol{u}_{h})$, both described in Section 3.3.1 below.

• •

Concerning $\left(q^{+}\,\widehat{c},z_{h}\right)_{0,\Omega}$, this term was approximated by

[TABLE]

• •

Regarding the mixed problem, the term $\mathcal{A}(c_{h};\boldsymbol{u}_{h},\boldsymbol{v}_{h})$ was substituted by

[TABLE]

with local forms

[TABLE]

where, similarly as before, $S^{K}_{\mathcal{A}}(\cdot,\cdot)$ is a stabilization term and $\nu_{\mathcal{A}}^{K}(c_{h})$ a constant, both described in Section 3.3.1 below.

• •

Finally, the term $(\boldsymbol{\gamma}(c_{h}),\boldsymbol{v}_{h})_{0,\Omega}$ was replaced by

[TABLE]

At this point, we highlight that the bilinear form $B(\cdot,\cdot)$ needs not to be substituted since it is computable for VE functions due to the choice of degrees of freedom (20). Furthermore, the right hand side term $\left(G,q_{h}\right)_{0,\Omega}$ remains unchanged.

Remark 4*.*

Note that we here use the convention that terms which are written in caligraphic letters, such as $\mathcal{M}_{h}$, $\mathcal{D}_{h}$ and $\mathcal{A}_{h}$, include a stabilization term, whereas those in non-caligraphic fashion and those of the form $(\cdot,\cdot)_{h}$ with subscript $h$ do not. In general, the terms of the type $(\cdot,\cdot)_{h}$ are approximations of the corresponding (possibly weighted) $L^{2}$ scalar products $(\cdot,\cdot)_{0,\Omega}$, obtained by introducing projections onto polynomials for all virtual functions, but not for the data terms that are known exactly.

3.3.1 Construction of the stabilizations

Here, we specify the assumptions on the stabilizations $S^{K}_{\mathcal{M}}(\cdot,\cdot):\,Z_{h}\times Z_{h}\to\mathbb{R}$, $S^{K}_{D}(\cdot,\cdot):\,Z_{h}\times Z_{h}\to\mathbb{R}$, and $S^{K}_{\mathcal{A}}(\cdot,\cdot):\,\boldsymbol{V}_{h}\times\boldsymbol{V}_{h}\to\mathbb{R}$, in (25), (28) and (30), respectively.

We require that these terms represent computable, symmetric, and positive definite bilinear forms that satisfy, for all $K\in\mathcal{T}_{h}$, the following property: there exist positive constants $M_{0}^{\mathcal{M}}$, $M_{1}^{\mathcal{M}}$, $M_{0}^{\mathcal{D}}$, $M_{1}^{\mathcal{D}}$, $M_{0}^{\mathcal{A}}$, $M_{1}^{\mathcal{A}}$, which are independent of $h$ and $K$, such that

[TABLE]

Note that continuity follows immediately from the properties:

[TABLE]

for all $z_{h},\widetilde{z_{h}}\in Z_{h}\cap\ker(\Pi^{0,K}_{k+1})$; analogously for the other forms. In practice, under mesh assumptions (D1)-(D2), one can take the following scaled stabilizations corresponding to the degrees of freedom:

[TABLE]

Regarding the constants appearing in front of the stabilizations in (25), (28) and (30), respectively, we pick:

[TABLE]

where $\Pi^{0,K}_{0}:\,L^{2}(K)\to\mathbb{P}_{0}(K)$ and $\boldsymbol{\Pi^{0,K}_{0}}:\,[L^{2}(K)]^{2}\to[\mathbb{P}_{0}(K)]^{2}$ are the $L^{2}$ projectors onto scalar and vector valued constants, respectively.

3.3.2 Well-posedness of the semidiscrete problem

We first define the constants

[TABLE]

Analogously, we introduce $\nu_{\mathcal{D}}^{-}$, $\nu_{\mathcal{D}}^{+}$, $\nu_{\mathcal{A}}^{-}$ and $\nu_{\mathcal{A}}^{+}$. Recalling (3) and (6), it is easy to check the following (mesh-uniform) bounds for the above constants:

[TABLE]

Then, similarly as for their continuous counterparts, the following continuity and coercivity properties for $\mathcal{M}_{h}(\cdot,\cdot)$ $\mathcal{D}_{h}(\cdot;\cdot,\cdot)$, and $\mathcal{A}_{h}(\cdot;\cdot,\cdot)$, defined in (25), (28) and (30), respectively, hold true.

Lemma 3.2**.**

For $\mathcal{M}_{h}(\cdot,\cdot)$, it holds, for all $c_{h},z_{h}\in Z_{h}$,

[TABLE]

Concerning $\mathcal{D}_{h}(\cdot;\cdot,\cdot)$, this form satisfies, for all $\boldsymbol{u}_{h}\in\boldsymbol{V}_{h}$ and $c_{h},z_{h}\in Z_{h}$,

[TABLE]

Regarding $\mathcal{A}_{h}(\cdot;\cdot,\cdot)$, for all $c_{h}\in Z_{h}$ and $\boldsymbol{u}_{h},\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}$, it yields

[TABLE]

Thus, $\mathcal{A}_{h}(c_{h};\cdot,\cdot)$ is coercive on the kernel

[TABLE]

with respect to $\lVert\cdot\rVert_{\boldsymbol{V}_{h}}$, where $\mathcal{K}$ is given in (13).

Proof.

The continuity bound in (35) follows directly by using

[TABLE]

and then estimating

[TABLE]

where the Pythagorean theorem was applied in the last equality. For the coercivity bound, one can use (6), (32), and the Pythagorean theorem.

Regarding the continuity estimate for $\mathcal{D}_{h}(\cdot;\cdot,\cdot)$, by using a splitting of the form (39), together with an estimate as in (14), one can deduce at the local level

[TABLE]

By application of a polynomial inverse estimate [26, Lemma 4.5.3], the continuity of the $L^{2}$ projector, and the Hölder inequality, we further estimate

[TABLE]

After inserting (41) into (40), taking the splitting into account, and summing over all elements, the stated bound follows. Concerning the coercivity bound for $\mathcal{D}_{h}(\cdot,\cdot)$, one can proceed similarly as in (16) for the consistency part, and employ (32) for the stabilization term, to obtain elementwise

[TABLE]

We now note that the definitions of $\Pi^{\nabla,K}_{k+1}$ and $\boldsymbol{\Pi^{0,K}_{k}}$ easily yield

[TABLE]

The estimate then follows with (42), the Pythagorean theorem and summation over all elements.

The estimates for $\mathcal{A}_{h}(\cdot;\cdot,\cdot)$ are derived in a similar fashion as those for $\mathcal{M}_{h}(\cdot,\cdot)$, using (6). The coercivity on $\mathcal{K}_{h}$ follows from the fact

[TABLE]

owing to the definition of $\boldsymbol{V}_{h}(K)$ in (18). ∎

Well-posedness of problem (24) can be shown by combining the results in [58] for parabolic problems with those in [27, 9] for mixed problems, using Lemma 3.2. More precisely, in the spirit of the two-step strategy applied in [37] for FEM, one can first show that for any given $c_{h}(t)\in L^{\infty}(\Omega)$, $t\in J$, the mixed problem

[TABLE]

admits a unique solution by applying the techniques in [27, 9], and then, by using the Gronwall lemma and Picard-Lindelöf (see e.g. [23, Ch.1.10]), that $c_{h}(t)$ is uniquely determined by the discrete concentration equation

[TABLE]

see also [58]. We do not write here the details since we focus directly on the fully discrete case, see the next section.

3.4 Fully discrete formulation

Here, our goal is to formulate a fully discrete version of (24).

To start with, we introduce a sequence of time steps $t_{n}=n\tau$, $n=0,\dots,N$, with time step size $\tau$. Next, we define $\boldsymbol{u}^{n}:=\boldsymbol{u}(t_{n})$, $p^{n}:=p(t_{n})$, $c^{n}:=c(t_{n})$, $G^{n}:=G(t_{n})$, $(q^{+})^{n}:=q^{+}(t_{n})$, and $\widehat{c}^{\,n}:=\widehat{c}(t_{n})$ as the evaluations of the corresponding functions at time $t_{n}$, $n=0,\dots,N$. Moreover, we denote by $\boldsymbol{U}^{n}\approx\boldsymbol{u}_{h}(t_{n})$, $P^{n}\approx p_{h}(t_{n})$ and $C^{n}\approx c_{h}(t_{n})$, the approximations of the semidiscrete solutions at those times when using a time integrator method. Among many time discretization schemes, we here make a computationally cheap choice by choosing a backward Euler method that is explicit in the nonlinear terms. The fully discrete system consequently reads as follows:

• •

for $n=0$: Given $c_{0,h}\in Z_{h}$, solve

[TABLE]

for all $\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}$ and $q_{h}\in Q_{h}$.

• •

for $n=1,\dots,N$: Solve first the concentration equation for $C^{n}$:

[TABLE]

for all $z_{h}\in Z_{h}$, where $C^{0}:=c_{0,h}$. Then, solve the mixed problem for $\boldsymbol{U}^{n}$ and $P^{n}$:

[TABLE]

for all $\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}$ and $q_{h}\in Q_{h}$.

Lemma 3.3**.**

Given $\tau>0$, provided that $G^{n},(q^{+})^{n},P^{n},C^{n}\in L^{\infty}(\Omega)$, $\boldsymbol{\gamma}(C^{n})\in[L^{2}(\Omega)]^{2}$, and $\boldsymbol{U}^{n}\in[L^{\infty}(\Omega)]^{2}$, for all $n=0,\dots,N$, the formulation (43)-(45) is uniquely solvable.

Proof.

Similarly as for the semidiscrete case, well-posedness of (43) and (45) follows by using the tools of [27, 9]. Regarding (44), we first rewrite that equation as

[TABLE]

We observe that all of the term are continuous with respect to the norm $\lVert\cdot\rVert_{1,\mathcal{T}_{h}}$. More precisely, for $\mathcal{M}_{h}(\cdot,\cdot)$ and $\mathcal{D}_{h}(\boldsymbol{U}^{n-1};\cdot,\cdot)$, continuity follows from Lemma 3.2 and the definition of the broken $H^{1}$ norm. Next, for the term involving $(q^{+})^{n}$, we simply apply the Cauchy-Schwarz inequality and the stability of the $L^{2}$ projector. Finally, for the term with $\Theta_{h}$, we estimate

[TABLE]

where we also employed an inverse inequality as in (41). Thus, by the Lax-Milgram lemma, it only remains to show that the left hand side of (46) is coercive with respect to $\lVert\cdot\rVert_{1,\mathcal{T}_{h}}$. This is however a direct consequence of

[TABLE]

owing to the fact that $q^{+}$ and $q^{-}$ are non-negative, and the coercivity bounds (35) and (36). ∎

Note that both problems (44) and (45) represent linear systems of equations which are decoupled from each other in the sense that, first, given $\widehat{c}^{\,n}$ and $(q^{+})^{n}$, one can determine $C^{n}$ with knowledge of $\boldsymbol{U}^{n-1}$ only, and then one can use $C^{n}$ to compute $\boldsymbol{U}^{n}$ and $P^{n}$. The quantity $P^{n}$ does in fact not influence the calculation of $C^{n}$ directly, but rather takes the role of a Lagrange multiplier and derived variable. This decoupling, combined with the fact that the systems to be solved at each time step are linear, makes the method quite cheap per iteration.

4 Error analysis for the fully discrete problem

The error analysis is performed in two steps: firstly, we estimate the discretization errors for the velocity and pressure, $\lVert\boldsymbol{u}^{n}-\boldsymbol{U}^{n}\rVert_{0,\Omega}$ and $\lVert p^{n}-P^{n}\rVert_{0,\Omega}$, respectively, and then, in the second step, the concentration error $\lVert c^{n}-C^{n}\rVert_{0,\Omega}$. In the following analysis, we assume all the needed regularity of the exact solution. Although such high regularity will not often be available in practice, the purpose of the following analysis is to give a theoretical backbone to the proposed scheme and to investigate its potential accuracy in the most favorable scenario.

4.1 An auxiliary result

The subsequent technical lemma will serve as an auxiliary result in the derivation of the error estimates and will be used in several occasions.

Lemma 4.1**.**

Let $r,s,t\in\mathbb{N}_{0}$. Denote by $\Pi^{0}_{r}$ and $\boldsymbol{\Pi^{0}_{s}}$, the elementwise defined $L^{2}$ projectors onto scalar and vector valued polynomials of degree at most $r$ and $s$, respectively. Given a scalar function $\sigma\in H^{m_{r}}(\mathcal{T}_{h})$, $0\leq m_{r}\leq r+1$, let $\kappa(\sigma)$ be a tensor valued piecewise Lipschitz continuous function with respect to $\sigma$. Further, let $\widehat{\sigma}\in L^{2}(\Omega)$, and let $\boldsymbol{\chi}$ and $\boldsymbol{\psi}$ be vector valued functions. We assume that $\kappa(\sigma)\in[L^{\infty}(\Omega)]^{2\times 2}$, $\boldsymbol{\chi}\in[H^{m_{s}}(\mathcal{T}_{h})\cap L^{\infty}(\Omega)]^{2}$, $\boldsymbol{\psi}\in[L^{2}(\Omega)]^{2}$, and $\kappa(\sigma)\boldsymbol{\chi}\in[H^{m_{t}}(\mathcal{T}_{h})]^{2}$, for some $0\leq m_{s}\leq s+1$ and $0\leq m_{t}\leq t+1$. Then,

[TABLE]

Proof.

We first write

[TABLE]

Then, for the first part on the right hand side of (47), we recall that $\boldsymbol{\Pi^{0}_{t}}$ is an $L^{2}$ projection and derive, on each element $K\in\mathcal{T}_{h}$,

[TABLE]

where in the last step we used Lemma 3.1 and the fact that $\kappa$ is Lipschitz continuous with respect to $\sigma$. The term $\lVert\boldsymbol{\Pi^{0,K}_{s}}\boldsymbol{\chi}\rVert_{\infty,K}$ is estimated as in (41). Concerning the second part on the right hand side of (47), we have, for each $K\in\mathcal{T}_{h}$,

[TABLE]

where we used again the Lipschitz continuity of $\kappa$, the continuity properties of the $L^{2}$ projectors, and the bound (41). The assertion of the lemma follows after combining the estimates and summing over all elements. ∎

Note that the above lemma can be easily transferred to the cases where $\sigma$, $\kappa(\sigma)$, $\chi$, and $\psi$ are scalar, and to vector valued $\boldsymbol{\sigma}$, $\boldsymbol{\chi}$ and scalar $\kappa(\boldsymbol{\sigma})$, $\psi$.

In the special case of $\chi=1$ and vector valued $\boldsymbol{\kappa}$, an adaptation of Lemma 4.1 gives

[TABLE]

4.2 Bounds on $\lVert\boldsymbol{u}^{n}-\boldsymbol{U}^{n}\rVert_{0,\Omega}$ and $\lVert p^{n}-P^{n}\rVert_{0,\Omega}$

We consider the mixed problem

[TABLE]

where $C^{n}\in Z_{h}$ is the numerical solution of the concentration equation (44) for $n=1,\dots,N$, and $C^{0}=c_{0,h}$. The goal is to derive an upper bound for $\lVert\boldsymbol{u}^{n}-\boldsymbol{U}^{n}\rVert_{0,\Omega}$ and $\lVert p^{n}-P^{n}\rVert_{0,\Omega}$ with respect to $\lVert c^{n}-C^{n}\rVert_{0,\Omega}$. For the analysis, we basically follow the ideas of [27, 9] with the major differences that, here, $\mathcal{A}_{h}(C^{n};\cdot,\cdot)$ is not consistent with respect to $[\mathbb{P}_{k}(K)]^{2}$ due to presence of $C^{n}$, and, additionally, the right hand side of (49) is inhomogeneous.

Theorem 4.2**.**

Given $C^{n}\in Z_{h}$, let $(\boldsymbol{U}^{n},P^{n})\in\boldsymbol{V}_{h}\times Q_{h}$ be the solution to (49). Let us assume that for the exact solution $(\boldsymbol{u}^{n},p^{n},c^{n})$ to (10) at time $t_{n}$, it holds $\boldsymbol{u}^{n}\in[H^{k+1}(\mathcal{T}_{h})]^{2}$, $p^{n}\in H^{k+1}(\mathcal{T}_{h})$, and $c^{n}\in H^{k+1}(\mathcal{T}_{h})$. Furthermore, we suppose that $\boldsymbol{\gamma}(c)$ and $A(c)$ are piecewise Lipschitz continuous functions with respect to $c\in L^{2}(\Omega)$, and that $\boldsymbol{\gamma}(c^{n}),A(c^{n})\boldsymbol{u}^{n}\in[H^{k+1}(\mathcal{T}_{h})]^{2}$. Then, the following error estimates hold for all $k\in\mathbb{N}_{0}$:

[TABLE]

where $\zeta_{1}^{n}$-$\zeta_{4}^{n}$ are positive constants independent of $h$ and depending only on the specified functions.

Proof.

The estimate for $\lVert\boldsymbol{U}^{n}-\boldsymbol{u}^{n}\rVert_{0,\Omega}$ can be obtained as follows.

By using the second equality in (49), we have $\operatorname{\rm div}\boldsymbol{U}^{n}=\Pi^{0}_{k}G^{n}$ (use that $\operatorname{\rm div}\boldsymbol{U}^{n}\in\mathbb{P}_{k}(K)$ for every $K\in\mathcal{T}_{h}$), where we recall that $(\Pi^{0}_{k})_{|_{K}}=\Pi^{0,K}_{k}$. Define now the interpolant $\boldsymbol{u}^{n}_{I}\in\boldsymbol{V}_{h}$ via the degrees of freedom (20):

[TABLE]

Then, it holds [9, eq.(28)]

[TABLE]

Moreover, one has $\operatorname{\rm div}\boldsymbol{u}^{n}_{I}=\Pi^{0}_{k}G^{n}$. Thus, setting $\boldsymbol{\delta}^{n}:=\boldsymbol{U}^{n}-\boldsymbol{u}^{n}_{I}$, it holds that $\boldsymbol{\delta}^{n}\in\mathcal{K}_{h}\subset\mathcal{K}$, where $\mathcal{K}_{h}$ and $\mathcal{K}$ were defined in (38) and (13), respectively, and therefore, $\lVert\boldsymbol{\delta}^{n}\rVert_{\boldsymbol{V}_{h}}=\lVert\boldsymbol{\delta}^{n}\rVert_{0,\Omega}$. Owing to the assumptions on $a(\cdot)$ in (6) together with (37), we have, further using (49) with $\boldsymbol{v}_{h}=\boldsymbol{\delta}^{n}\in\mathcal{K}_{h}$ and (10),

[TABLE]

The terms $T_{1}$-$T_{3}$ are bounded as follows:

• •

term $T_{1}$: We use equation (48) with $\boldsymbol{\kappa}=\boldsymbol{\gamma}$, $\sigma=c^{n}$, $\widehat{\sigma}=C^{n}$, $\boldsymbol{\psi}=\boldsymbol{\delta}^{n}$, $r=k+1$, $t=k$, and $m_{r}=m_{t}=k+1$, and obtain

[TABLE]

• •

term $T_{2}$: Owing to the continuity properties (37) of $\mathcal{A}_{h}(\cdot;\cdot,\cdot)$ and the interpolation error estimate (50), it holds

[TABLE]

• •

term $T_{3}$: We have

[TABLE]

For the term $T_{3}^{A}$, we use Lemma 4.1 with $\kappa=A$, $\sigma=c^{n}$, $\widehat{\sigma}=C^{n}$, $\boldsymbol{\chi}=\boldsymbol{u}^{n}$, $\boldsymbol{\psi}=\boldsymbol{\delta}^{n}$, $r=k+1$, $s=t=k$, and $m_{r}=m_{s}=m_{t}=k+1$, to get

[TABLE]

On the other hand, the term $T_{3}^{B}$ can be bounded with (32), (6), and Lemma (3.1):

[TABLE]

After plugging the bounds obtained for $T_{1}$-$T_{3}$ into (51), dividing by $\lVert\boldsymbol{\delta}^{n}\rVert_{0,\Omega}$, using the triangle inequality in the form

[TABLE]

and employing (50), the convergence result follows.

The error estimate for the term $\lVert P^{n}-p^{n}\rVert_{0,\Omega}$ follows easily by combining the above ideas with the argument in [27, Theorem 6.1] and is therefore not shown. ∎

4.3 Bounds on $\lVert c^{n}-C^{n}\rVert_{0,\Omega}$

For fixed $\boldsymbol{u}(t)\in\boldsymbol{V}$ and $t\in J$, we define the projector $\mathcal{P}_{c}:\,Z\to Z_{h}$ (that to each $c\in Z$ associates $\mathcal{P}_{c}c\in Z_{h}$) by

[TABLE]

for all $z_{h}\in Z_{h}$, where

[TABLE]

with

[TABLE]

Lemma 4.3**.**

The projector $\mathcal{P}_{c}:\,Z\to Z_{h}$ given in (52) is well-defined under the assumption that $\boldsymbol{u}$, $q^{+}$, and $q^{-}$ are bounded in $L^{\infty}(\Omega)$ for all $t\in J$.

Proof.

By the Lax-Milgram lemma, we have to show that the left hand side of (52) defines a continuous and coercive bilinear form and that the right hand side is a continuous functional with respect to $\lVert\cdot\rVert_{1,\mathcal{T}_{h}}$. Continuity of the latter one is obtained by combining (14) with

[TABLE]

By using (36) and performing similar computations as in the proof of Lemma 3.3, continuity of $\Gamma_{c,h}$ follows:

[TABLE]

where $\zeta$ only depends on the specified functions. Regarding the coercivity of $\Gamma_{c,h}$, we first estimate

[TABLE]

where we recall that $(\Pi^{0}_{k+1})_{|_{K}}=\Pi^{0,K}_{k+1}$ for all $K\in\mathcal{T}_{h}$. Then, combining this result with (36) yields

[TABLE]

with $\overline{z_{h}}$ denoting the $L^{2}(\Omega)$ projection of $z_{h}$ onto $\mathbb{P}_{0}(\Omega)$. Since $\overline{z_{h}}$ coincides with the average of $z_{h}$, one can use a Poincaré-Friedrichs inequality, see e.g. [24], to deduce

[TABLE]

and consequently the coercivity of $\Gamma_{c,h}$. ∎

Lemma 4.4**.**

We assume that $\boldsymbol{u}\in[H^{k+1}(\mathcal{T}_{h})\cap L^{\infty}(\Omega)]^{2}$, $c\in H^{k+2}(\mathcal{T}_{h})\cap W^{1,\infty}(\mathcal{T}_{h})$, $q^{+},q^{-}\in L^{\infty}(\Omega)$, $(q^{+}+q^{-})c\in H^{k+1}(\mathcal{T}_{h})$, $\boldsymbol{u}\,c\in[H^{k+1}(\mathcal{T}_{h})]^{2}$, $\boldsymbol{u}\cdot\nabla c\in H^{k+1}(\mathcal{T}_{h})$, and $D(\boldsymbol{u})\nabla c\in[H^{k+1}(\mathcal{T}_{h})]^{2}$ for all $t\in J$. Then, the following error bounds for $c-\mathcal{P}_{c}c$, where $\mathcal{P}_{c}c$ is defined in (52), hold for all $k\in\mathbb{N}_{0}$:

[TABLE]

where the constants $\xi_{1},\xi_{0}>0$ only depend on the listed terms and are independent of $h$.

Proof.

We focus on the error estimate in the broken $H^{1}$ norm at a fixed time $t\in J$. First, we state the following result. Given $c\in H^{k+2}(\mathcal{T}_{h})$, there exists an interpolant $c_{I}\in Z_{h}$ such that the following bounds hold true (see for instance [30, 15, 25]):

[TABLE]

After denoting $\nu:=\mathcal{P}_{c}c-c_{I}$, one obtains with the coercivity of $\Gamma_{c,h}$, see the proof of Lemma 4.3, and the definition of $\mathcal{P}_{c}c$ in (52),

[TABLE]

for a constant $M>0$. By employing the definitions of $\Gamma_{c}$ and $\Gamma_{c,h}$ in (53), the term $S_{1}$ is split as follows:

[TABLE]

For $S_{1}^{A}$, we have

[TABLE]

where in the inequality we applied Lemma 4.1 to estimate the first part on the right hand side of $S_{1}^{A}$, and made use of the continuity properties (32) of $S^{K}_{D}(\cdot,\cdot)$, the trivial continuity property of $\Pi^{\nabla,K}_{k+1}$ in the $H^{1}$ seminorm and its approximation properties (stated in Lemma 3.1) to estimate the stabilization term.

Next, for $S_{1}^{B}$, we compute

[TABLE]

where in the last inequality we used Lemma 4.1 with $\kappa=id$ and $\sigma=\boldsymbol{u}$ for the first and third term inside the curly bracket, and $\kappa=q^{+}+q^{-}$ and $\sigma=1$ for the second one.

Finally, for $S_{1}^{C}$, it holds with the definition of the $L^{2}$ projector and Lemma 3.1

[TABLE]

On the other hand, for $S_{2}$, we use the continuity of $\Gamma_{c,h}$ in (54), together with the interpolation error estimate (56), to derive

[TABLE]

The error bound in the broken $H^{1}$ norm follows by plugging first the estimates for $S_{1}^{A}$, $S_{1}^{B}$, and $S_{1}^{C}$ into $S_{1}$, then those obtained for $S_{1}$ and $S_{2}$ into (57), using the definition of the $H^{1}$ norm, dividing by $\lVert\nu\rVert_{1,\mathcal{T}_{h}}$, and using the triangle inequality in the form

[TABLE]

together with the approximation properties (56) of the interpolant $c_{I}$.

The $L^{2}$ error bound can be derived by combining the above arguments with a standard duality argument as in [58], also recalling the convexity of $\Omega$; it is omitted here. ∎

By differentiation of (52) in time and use of similar techniques as in the proof of Lemma 4.4, an analogous result can be obtained for $\frac{\partial}{\partial t}(c-\mathcal{P}_{c}c)$, summarized in the following corollary.

Corollary 4.5**.**

Provided that the continuous data and solution are sufficiently regular in space and time, it holds

[TABLE]

where the constants $\xi_{1,t},\xi_{0,t}>0$ are independent of $h$.

Moreover, we will later on need the two subsequent bounds.

Lemma 4.6**.**

Under sufficient smoothness of the continuous data and solution, it holds

[TABLE]

where $\xi_{0,t}$ can be found in Corollary 4.5.

Proof.

We estimate

[TABLE]

The term $(I)$ can be estimated exactly as for standard finite elements, see for instance [56]:

[TABLE]

where we also applied the Hölder inequality in the last step. Concerning $(II)$, this term can be bounded as follows, using Corollary 4.5:

[TABLE]

The statement of the lemma follows. ∎

Lemma 4.7**.**

Provided that the continuous data and solution are sufficiently regular in space and time, it holds

[TABLE]

where $\zeta_{1}^{n-1}$ and $\zeta_{2}^{n-1}$ are the constants from Theorem 4.2.

Proof.

By using the triangle inequality, one obtains

[TABLE]

The first term on the right hand side is estimated by

[TABLE]

and the second one term is bounded with Theorem 4.2. ∎

Now, we have all the ingredients to bound $\lVert c^{n}-C^{n}\rVert_{0,\Omega}$.

Theorem 4.8**.**

Let the mesh assumptions (D1)-(D3) be satisfied. Then, provided that the continuous data and solutions are sufficiently regular, it yields

[TABLE]

where the regularity terms $\varphi_{1},\varphi_{2}$ and the positive constant $\eta$ now depend on $\boldsymbol{u}$, $c$, $q^{+}$, $q^{-}$, $\widehat{c}$, $\frac{\partial\boldsymbol{u}}{\partial t}$, $\frac{\partial^{2}\boldsymbol{u}}{\partial t^{2}}$, $\frac{\partial c}{\partial t}$, and $\frac{\partial^{2}c}{\partial t^{2}}$ (and products of these functions).

Proof.

To start with, we write

[TABLE]

Equation (55) gives a bound on $\rho^{n}$. In order to deal with $\vartheta^{n}$, we use the continuous concentration equation (12) with $z=\vartheta^{n}$, the fully discretized version (44) with $z_{h}=\vartheta^{n}$, and the definition of the projector $\mathcal{P}_{c}c^{n}$ in (52) with $z_{h}=\vartheta^{n}$:

[TABLE]

Owing to the coercivity properties in (36), the second term on the left hand side of (58) can be estimated by

[TABLE]

with some constant $D_{\ast}>0$ independent of $h$ and $\boldsymbol{U}^{n-1}$.

The terms $R_{1}$-$R_{5}$ on the right hand side of (58) are estimated as follows:

• •

term $R_{1}$: Using the definition of $\mathcal{M}_{h}(\cdot,\cdot)$ in (25) together with (32) yields

[TABLE]

The term $R_{1}^{A}$ is estimated by using the continuity of the $L^{2}$ projector, the assumption (6) on $\phi$, and the approximation properties in Lemma (3.1):

[TABLE]

Next, we bound $R_{1}^{B}$ with similar tools as for $R_{1}^{A}$:

[TABLE]

Thus, we deduce with Lemma 4.6

[TABLE]

with the obvious definitions for the regularity terms $R_{1}^{n,1}$, $R_{1}^{n,2}$, and $R_{1}^{n,3}$.

• •

term $R_{2}$: By the definition of $\Theta_{h}(\cdot;\cdot,\cdot)$ in (27), the identity $\vartheta^{n}=C^{n}-\mathcal{P}_{c}c^{n}$, and the fact that $(q^{+})^{n}$ and $(q^{-})^{n}$ are non-negative, it holds

[TABLE]

The above equation, after adding zero in the form

[TABLE]

to the right hand side, can be equivalently expressed as

[TABLE]

For $R_{2}^{A}$, we estimate

[TABLE]

We now use an inverse estimate [26, Lemma 4.5.3], the continuity of $\boldsymbol{\Pi^{0,K}_{k}}$, a triangle inequality, the assumption that $\mathcal{T}_{h}$ is quasi-regular, and Lemma 4.4, to deduce, for every $K\in\mathcal{T}_{h}$,

[TABLE]

Recalling Lemma 4.7, the definitions of $\vartheta^{n-1}$ and $\rho^{n-1}$, and Lemma 4.4, we get

[TABLE]

thus implying

[TABLE]

with the obvious definitions for the regularity terms $R_{2}^{n,1}$, $R_{2}^{n,2}$, and $R_{2}^{n,3}$.

The term $R_{2}^{B}$ can be bounded analogously to $R_{2}^{A}$, giving

[TABLE]

Using again the bound (62), one obtains

[TABLE]

Thus,

[TABLE]

• •

term $R_{3}$: We use the definition of $\mathcal{D}_{h}(\cdot;\cdot,\cdot)$ in (28), a standard Hölder inequality in the spirit of (15), the estimate (61), the scaling properties of the stabilization in (32), and the Lipschitz continuity of $D(\cdot;\cdot,\cdot)$ and $\nu_{D}^{K}$ in (34), to deduce

[TABLE]

Hence, with (62) we have

[TABLE]

• •

term $R_{4}$: The use of Lemma 4.4 yields

[TABLE]

with the obvious definition of $R_{4}^{n,1}$.

• •

term $R_{5}$: The approximation properties in Lemma 3.1 yield

[TABLE]

with the obvious definition of $R_{5}^{n,1}$.

We now insert (59) and the bounds on $R_{1}$-$R_{5}$ into (58). Afterwards, we observe that all regularity terms $\{R_{J}^{n,i}\}$ above only depend on the continuous solution and can be assumed to be bounded uniformly in $h$. We only keep track of the terms $R_{1}^{n,2}$ and $R_{1}^{n,3}$. This yields

[TABLE]

with the positive scalars

[TABLE]

Next, we introduce, for all $w_{h}\in Z_{h}$, the discrete norm

[TABLE]

Owing to Lemma 3.2, there exist positive constants $c_{\ast}$ and $c^{\ast}$, such that, for all $w_{h}\in Z_{h}$, it holds

[TABLE]

Reshaping (63), and employing (65) and (66), then gives

[TABLE]

The terms $T_{1}$ and $T_{2}$ are bounded as follows:

[TABLE]

where we used (39) and (65) in the first step. The term $T_{3}$ is bounded as follows

[TABLE]

Next, we plug (68) and (69) into (67), cancel the terms $\tau D_{\ast}|\vartheta^{n}|^{2}_{1,\mathcal{T}_{h}}$ and manipulate the resulting inequality, to obtain

[TABLE]

Moreover, we estimate

[TABLE]

Hence,

[TABLE]

Defining

[TABLE]

and solving the recursion then leads to

[TABLE]

where we recall that $n\leq T/\tau$ with $T$ the final time instant. With (66) the estimate in the $L^{2}$ norm is a direct consequence:

[TABLE]

The initial term $\lVert\vartheta^{0}\rVert_{0,\Omega}^{2}$ is estimated by

[TABLE]

where we applied Lemma 4.4. Moreover, using (64), the fact that $\sum_{j=1}^{n}\tau\leq T$, and the definitions of $R_{1}^{j,2}$ and $R_{1}^{j,3}$ in (60), after some simple manipulations, we obtain

[TABLE]

The assertion of the theorem follows by combining (70) with (72) and (71). ∎

5 Numerical experiments

In this section, we demonstrate the performance of the method on the basis of numerical experiments, focusing on the lowest order case $k=0$. To this purpose, we first consider an ideal test case (Example 1), and then a more realistic one (Example 2). The aim of the first test is to validate (also numerically) the convergence of the method on a problem with regular known solution, whereas those of the second test is to check the method’s performance on a well-known benchmark that mimics a more realistic situation.

Example 1: Here, we study a generalized version of (1), given by

[TABLE]

endowed with the boundary and initial conditions in (4) and (5), respectively. We fix $\Omega=(0,1)^{2}$ and pick the same choice of parameters as in [44], namely $T=0.01$, $\phi=1$, $D(\boldsymbol{u})=|\boldsymbol{u}|+0.02$, $d_{m}=0.02$, $d_{\ell}=d_{t}=1$, $c_{0}=0$, $\boldsymbol{\gamma}(c)=0$, and $a(c)=(c+2)^{-1}$, where $f$ and $g$ are taken in accordance with the analytical solutions

[TABLE]

Plots of the exact solution at the final time $T$ are shown in Figures 1 and 2.

We employ a sequence of regular Cartesian meshes and Voronoi meshes, as portrayed in Figure 3. In addition to the current version, we also test the method when replacing the stabilization terms in (26), (29), and (31) by alternative ones:

[TABLE]

The alternative (diagonal) stabilizations are given by

[TABLE]

with

[TABLE]

where $\{\varphi_{j}^{K}\}_{\ell=1}^{\textrm{dim}Z_{h}(K)}$ and $\{\boldsymbol{\psi}_{j}^{K}\}_{\ell=1}^{{\textrm{dim}{V}_{h}(K)}}$ denote the local canonical basis functions for $Z_{h}(K)$ and $\boldsymbol{V}_{h}(K)$, and $\sigma>0$ is a safety parameter. In the forthcoming experiments, we set $\sigma=1e-3$. We highlight that these stabilizations are in fact modifications of the so-called D-recipe, which was introduced in [14] and has already been successfully applied in some variants to other model problems, such as the Helmholtz problem [50]. The first entry inside the max is simply the “diagonal part” of the consistency term of the local approximate forms in (26), (29), and (31), respectively, whereas the second terms correspond to the original stabilizations associated to the degrees of freedom in (33) multiplied by $\sigma$, which acts as a positivity safeguard. Importantly, it is easy to check that the error analysis can be easily extended to the new choice of stabilizations.

Due to the virtuality of the basis functions, we measure the following relative $L^{2}$ errors:

[TABLE]

where $C^{n}$, $\boldsymbol{U}^{n}$, and $P^{n}$ are the numerical solutions at the final time $T$.

The relative $L^{2}$ discretization errors for the concentration are plotted in Figure 4 in terms of the mesh size $h$ for both families of meshes and both variants of stabilizations. In order to better underline the expected linear convergence of the method both in $h$ and $\tau$ (see Theorem 4.8, recalling that $k=0$), the time step $\tau$ is chosen proportional to $h$. In other words, starting with the coarsest mesh and $\tau=T/5$, each subsequent case is obtained by dividing both $h$ (adopting a finer mesh) and $\tau$ by a factor of 2. Analogous plots are shown for the velocity and pressure variable errors in Figures 5 and 6. In all cases, the linear convergence rates are in accordance with Theorem 4.2 and Theorem 4.8. For the pressure discretization error, since the initial meshes are very coarse, we observe some pre-asymptotic regime when employing the original stabilizations in (33). This effect, however, is not present for the alternative stabilizations in (74). Both variants lead to similar results for the concentration and velocity errors.

Since the concentration often evolves more rapidly than the velocity and pressure, it could be worth to consider a cheaper variant of the discrete scheme (43)-(45) where the discrete velocity-pressure pair is updated only every R time steps (with $R\in{\mathbb{N}}$). This leads to a smaller number of linear system resolutions (possibly with a small reduction in accuracy) since only system (44) is solved at every time step, while (45) is solved only every R steps. In order to test this, we tried to run the same test above and compare the original version with the cheaper version with $R=5$. The difference in error was only at the fourth meaningful digit; we do not plot the graphs since these would completely overlap the ones of the original method.

Example 2: Next, we investigate the behavior of the method for Test 1 and Test 2 in [59, 33].

The problem is given in the form (1) with boundary conditions (4) and initial condition (5) over the spatial domain $\Omega=(0,1000)^{2}$ ft2. Moreover, $T=3600$ days and $\tau=36$ days. At the upper right corner, i.e. at $[1000,1000]$, fluid with concentration $\widehat{c}=1.0$ is injected with rate $q^{+}=30$ ft2/day, whereas at the lower left corner, i.e. at $[0,0]$, material is absorbed with rate $q^{-}=30$ ft2/day. Both wells are henceforth treated as Dirac masses, which is admissible at the discrete level since the discrete functions are piecewise regular (which can be interpreted as an approximation of the Dirac delta by a localized function with support within the corner element and unitary integral). Furthermore, the following choices for the parameters are picked: $\phi=0.1$, $d_{\ell}=50$, $d_{t}=5$, $c_{0}=0$, $\boldsymbol{\gamma}(c)=0$, and $a(c)=80(1+(M^{\frac{1}{4}}-1)c)^{4}$, where

[TABLE]

Whereas $a(c)$ is constant for Test A, it changes rapidly across the fluid interface for Test B (which is in fact not covered by the theoretical analysis since $d_{m}=0$, but is interesting to study numerically) resulting in a much faster propagation of the fluid concentration front along the diagonal direction ($d_{\ell}\gg d_{t}$). This effect is known as macroscopic fingering phenomenon[39].

For this example, we used a regular 25x25 Cartesian mesh and we employed the more sophisticated stabilization in (74). Since Test B is highly convection-dominated, pure application of our method leads to local disturbances in the form of overshoots and undershoots of the numerical solution for the concentration, typical in the context of convection-dominated problem. To this purpose, for this test case, we employ the flux-corrected transport (FCT) algorithm with linearization [48, 47]. The FCT scheme with linearization for convection-dominated flow problems operates in two steps: (1) advance the solution in time by a low-order overly diffusive scheme to suppress spurious oscillations, (2) correct the solution using (linear) antidiffusive fluxes. In that way the computed solution does not show spurious oscillations and layers are not smeared.

Due to the fact that no analytical solutions are available for Test A and Test B, we plot the numerical solutions (and the corresponding contour plots) for the concentration after 3 and 10 years. These times correspond to $n=30$ and $n=100$, respectively. For visualization of the results, since the numerical solution is virtual but the nodal values are known, we simply add, inside each square, the barycenter with associated mean value of the nodal values, then create a triangulation based upon these points, and finally interpolate the function values linearly inside each triangle. In Figures 7 and 8, the results for Test A are portrayed, and in Figures 9 and 10, those for Test B. The results are similar to those obtained in [59, 33].

Acknowledgements

The first (L.B.d.V.) and last (G.V.) authors where partially supported by the European Research Council through the H2020 Consolidator Grant (grant no. 681162) CAVE, Challenges and Advancements in Virtual Elements. This support is gratefully acknowledged. The second author (A.P.) has been funded by the Austrian Science Fund (FWF) through the project P 29197-N32, and by the Doctoral Program (DK) through the FWF Project W1245.