# A conservative network element method for diffusion-advection-reaction problems

Julien Coatléven

► To cite this version:

Julien Coatléven. A conservative network element method for diffusion-advection-reaction problems. ESAIM: Mathematical Modelling and Numerical Analysis, 2023, 57 (4), pp.2007-2040. 10.1051/m2an/2023040. hal-04179433

HAL Id: hal-04179433

https://ifp.hal.science/hal-04179433

Submitted on 9 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A CONSERVATIVE NETWORK ELEMENT METHOD FOR DIFFUSION-ADVECTION-REACTION PROBLEMS

JULIEN COATLÉVEN*

**Abstract.** We derive a conservative network element method for heterogeneous and anisotropic diffusion problems by modifying the non-conservative version, and extend it to the approximation of an additional advection term. The numerical scheme possesses the flux formulation reminiscent of classical finite volume methods. Its convergence is naturally governed by the network element theory. Numerical results illustrate the good behavior of the method even on distorted point clouds.

**Mathematics Subject Classification.** 65N30, 65N12, 65N15.

Received January 13, 2023. Accepted April 28, 2023.

## INTRODUCTION

The huge literature on finite volume methods emphasizes their attractiveness in many scientific fields such as biology, plasma physics, hydrodynamics, oil reservoir simulation, etc.. They allow to deal easily with the non-linearity occurring in complex, coupled physical models while preserving relevant physical properties such as local mass conservation. For a quite recent review, we refer the reader to [18]. Multiple point flux approximations (MPFA, see [1–3]) are the natural generalization of the original two-point flux approximation (TPFA). They allow to construct consistent and cell-centered flux approximations on distorted meshes and lead to compact stencils, but suffer from conditional coercivity (see [4, 5]). New finite volume methods using additional face unknowns (the hybrid finite volume scheme, see [19, 21, 22]) or vertex unknowns (the vertex approximate gradient schemes VAG, see [8, 23–25]) have been introduced. They allow this time to obtain both coercivity and consistency on quite general, possibly nonconforming meshes, for strongly heterogeneous and anisotropic tensors. Finally, let us mention the virtual volume (or conservative first order virtual element) method of [10], that was derived in the virtual element [6, 15] framework and still lead to a finite volume method, that can be formally identified as a VAG scheme.

All those approaches are mesh-based, conservative meshless methods being much less common (see [7, 9, 28] for a review of meshless methods). The key point for conservation properties is to recover Stokes formula at the discrete level, also often named compatibility of the discretization. Among meshfree methods that allow it, we mention the uncertain grid method [17], the meshfree framework for conservation laws [37], the meshless volume scheme [26] and the meshfree mimetic method [34] which is a conservative version of the staggered meshless method of [35]. In a recent paper [13], the principles of a compatible meshfree method called the network element

*Keywords and phrases.* Meshless methods, variational method, network element method.

IFP Énergies nouvelles, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison, France.

*Corresponding author: julien.coatleven@ifpen.fr
method (NEM) were developed for heterogeneous and anisotropic diffusion problems. Convergence proofs and estimates were then obtained in a fairly general framework in [11, 13], while efficient ways of generating the scheme are presented in [14]. However, in its original form while compatible the NEM does not ensure local mass conservation as a finite volume scheme would, nor does it present an obvious flux-based formulation which would simplify its extension to some important non-linear flow problems. The objective of the present paper is two-fold: first we explain how to modify the non-conservative NEM method for heterogeneous and anisotropic diffusion-reaction problems of [13] to obtain a conservative version that retains the usual flux formulation of finite volume methods. Using the fact that the NEM is derived by mimicking virtual element principles, we will simply follow the algebraic approach of [10] that has already proved to be successful in the virtual element context. Secondly, we explain how to approximate advective terms with the NEM, recovering a NEM equivalent of the usual upwind scheme for advection familiar to finite volume practitioners.

Apart from allowing simpler discretization of domains with complex geometries by substituting point clouds to meshes, one of the main interest of meshless methods based on point networks is to ease the use of Lagrangian formulations of advection. Contrary to Lagrangian meshes where the quality of the deformed mesh is unlikely to be maintained over long periods of time, moving points can be done at almost zero cost without too much endangering solution quality. As a practical application of Lagrangian simulations where the NEM can be an interesting alternative to mesh-based methods, we mention as an illustrative example the so-called basin simulation problem, which is a perhaps less well known use of two phase porous media flow models than reservoir simulation but still an important one for oil and gas exploration. The idea of such models is to assess the hydrocarbon potential of a given geological prospect through the simulation of the entire petroleum system over millions of years (see [27] for details). Over such extended periods of time, sedimentary layers will be deposited, eroded, compacted, deformed and fractured. Tectonic uplift, compression, rifting and faulting will in particular produce complex movements that are hard to follow in an Eulerian setting. This is the reason why the sedimentary evolution is in general simulated in a Lagrangian framework, the fluid flow being then simulated in the moving sediment frame. This poses great meshing challenges: in particular, to avoid too much numerical diffusion and obtain a correct estimation of hydrocarbon quantities and position remeshing is in general avoided as much as possible. This leads to unique meshing challenges, as one should not only produce a mesh valid for a given complex faulted geometry at a given time, but one that will remain numerically valid while undergoing complex deformations for long periods of time. This remains up to this day a difficult, computationally expensive and not fully resolved issue [20, 29, 33, 36]. Meshless methods often here a potential alternative: by representing the geometry through point networks, they can undergo complex deformations while remaining a valid discretization.

The paper will be organized as follows: in a first section, we briefly recall the notion of discretization network and the associated approximate network geometry. Then, we explain how to derive a conservative network element approximation of a model advection-diffusion-reaction problem thanks to a discrete variational formulation. In a third section, we make some comments about convergence results and error estimates and in particular we explain that they can be obtained by a straightforward adaptation of their counterparts for the non conservative case. The fourth section is dedicated to deriving a flux formulation of the conservative NEM, that is easier to implement and makes more obvious the relationship with classical finite volumes on meshes. In the fifth section we propose a projection operator between networks, thus enabling resampling operations during simulation. Finally, in the last section we exhibit some numerical results to illustrate the good behavior of the method. We start by considering classical test cases for advection-diffusion-reaction of the literature, to explore the behavior of the NEM in the advection dominated and diffusion dominated regimes. Next, we consider a time-dependent advection diffusion problem on a moving domain, for which we describe a NEM scheme in both the Eulerian and Lagrangian frameworks. We conclude our numerical section by an application of the NEM to a simplified yet realistic basin simulation relying on an elementary two phase porous media flow model. Having a flux based NEM formulation the extension to this simple Darcy model is indeed quite straightforward.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. An example of network for a circular domain with two holes: blue circles are cells, orange triangle are interfaces, lines represent the network connectivity.</figcaption></figure>

# 1. DISCRETIZATION NETWORKS AND NETWORK GEOMETRIES

In this section, we recall for the reader's convenience the notions of discretization network and associated network geometry that were originally introduced in [12]. Let $\Omega$ be an open bounded connected subset of $\mathbb{R}^d$, $d \in \mathbb{N}\{0\}$, assumed to be at least Lipschitz. For any $\mathbf{x} \in \mathbb{R}^d$ and any $r > 0$, we denote $B(\mathbf{x}, r)$ the ball of radius $r$ centered at $\mathbf{x}$.

A discretization network $\mathcal{N}$ of $\Omega$ is defined from two sets of points $\mathcal{P}_T$ and $\mathcal{P}_F$, by setting $\mathcal{N} = \{\mathcal{T}, \mathcal{F}\}$ where the set of cells $\mathcal{T}$ is a set of pairs $K = \{\mathbf{x}_K, r_K\}$, with $\mathbf{x}_K \in \mathcal{P}_T$ strictly inside $\Omega$ and $r_K$ a strictly positive real number for any $K \in \mathcal{T}$ (see Fig. 1). We denote $h_K = 2r_K$. The set of interfaces, denoted $\mathcal{F}$, is a set of pairs $\sigma = \{\mathbf{x}_\sigma, T_\sigma\}$ with $\mathbf{x}_\sigma \in \mathcal{P}_F$ and $T_\sigma$ a subset of $\mathcal{T}$. It is subdivided into two subsets, the set of boundary interfaces $\mathcal{F}_{\text{ext}}$ and the set of interior interfaces $\mathcal{F}_{\text{int}}$. The set of boundary interfaces $\mathcal{F}_{\text{ext}}$ is such that for all $\sigma \in \mathcal{F}_{\text{ext}}$, $\mathbf{x}_\sigma$ is a point in $\cap_{K \in \mathcal{T}_\sigma} B(\mathbf{x}_K, r_K) \cap \partial\Omega$. The set of interior interfaces $\mathcal{F}_{\text{int}}$ is such that for all $\sigma \in \mathcal{F}_{\text{int}}$, $\mathbf{x}_\sigma$ is a point in $\cap_{K \in \mathcal{T}_\sigma} B(\mathbf{x}_K, r_K) \cap \bar{\Omega}$. We assume that $\Omega$ is correctly covered by the balls surrounding each cell point, i.e. $\bar{\Omega} \subset \cup_{K \in \mathcal{T}} B(\mathbf{x}_K, r_K)$, and also that for all $(K_1, K_2) \in \mathcal{T}^2$ such that $K_1 \neq K_2$, $\mathbf{x}_{K_1} \neq \mathbf{x}_{K_2}$. Finally, for any $K \in \mathcal{T}$ such that $\partial\Omega \cap \bar{B}(\mathbf{x}_K, r_K) \neq \emptyset$ there exists $\sigma \in \mathcal{F}_{\text{ext}}$ such that $K \in T_\sigma$, and for any $(K, L) \in \mathcal{T}^2$ such that $B(\mathbf{x}_K, r_K) \cap B(\mathbf{x}_L, r_L) \neq \emptyset$ there exists $\sigma \in \mathcal{F}$ such that $(K, L) \subset T_\sigma$. For any $K \in \mathcal{T}$, we denote $\mathcal{F}_K = \{\sigma \in \mathcal{F} \mid K \in T_\sigma\}$ the interfaces of $K$ implying that for any $\sigma \in \mathcal{F}$, $T_\sigma = \{K \in \mathcal{T} \mid \sigma \in \mathcal{F}_K\}$. We denote $B_K = B(\mathbf{x}_K, r_K)$ for any $K \in \mathcal{T}$, $h = \max_{K \in \mathcal{T}} h_K$ and $\mathbb{P}_k(\mathbb{R}^d)$ the set of polynomials of order $k$. A network is
said to be admissible if for any cell $K \in \mathcal{T}$, the set $(\boldsymbol{x}_{\sigma})_{\sigma \in \mathcal{F}_K}$ is unisolvent for first order polynomials (see [12] for details).

A network geometry is defined as a set of coefficients $\mathcal{G} = ((m_K)_{K \in \mathcal{T}}, (\eta_{K,\sigma})_{K \in \mathcal{T}, \sigma \in \mathcal{F}_K}, \varepsilon(\mathcal{G}))$ with

$$ \varepsilon(\mathcal{G}) = \left( (\varepsilon_{K}^{0,i})_{K \in \mathcal{T}, 1 \le i \le d}, (\varepsilon_{K}^{1,ij})_{K \in \mathcal{T}, 1 \le i,j \le d}, (\varepsilon_{\sigma}^{i})_{\sigma \in \mathcal{F}_{\text{int}}, 1 \le i \le d} \right). $$

The discrete measures $(m_K)_{K \in \mathcal{T}}$ are said to be admissible if and only if they satisfy

$$ m_K > 0 \quad \text{for all } K \in \mathcal{T}, \quad \text{and} \quad \sum_{K \in \mathcal{T}} m_K = |\Omega|, \qquad (1) $$

while the approximate consistency properties are given by

$$ \sum_{\sigma \in \mathcal{F}_K} \eta_{K, \sigma}^{i} = m_K \varepsilon_K^{0, i} \quad \forall K \in \mathcal{T}, \forall 1 \le i \le d, \qquad (2) $$

and

$$ \sum_{\sigma \in \mathcal{F}_K} \eta_{K, \sigma}^{i} (x_{\sigma}^{j} - x_{K}^{j}) = m_K (\delta_{ij} + \varepsilon_{K}^{1,ij}) \quad \forall K \in \mathcal{T}, \forall 1 \le i, j \le d, \qquad (3) $$

where $\delta_{ij}$ is the Kronecker symbol, and the approximate compatibility (or conservation) properties by

$$ \sum_{K \in \mathcal{T}_{\sigma}} \eta_{K \sigma}^{i} = \varepsilon_{\sigma}^{i} \quad \forall \sigma \in \mathcal{F}_{\text{int}}, \forall 1 \le i \le d. \qquad (4) $$

To measure the geometric approximation error $\varepsilon(\mathcal{G})$, we introduce the constants $\theta_A > 0$ and $p \ge 1$, both independent on $h$ and such that:

$$ |\varepsilon_K^{0,i}| \le \theta_A h_K^p \quad \text{and} \quad |\varepsilon_K^{1,ij}| \le \theta_A h_K^p \quad \forall K \in \mathcal{T}, \forall 1 \le i, j \le d, \qquad (5) $$

and

$$ |\varepsilon_{\sigma}^{i}| \le \theta_{A} \min_{K \in \mathcal{T}_{\sigma}} m_{K} h_{K}^{p} \quad \forall \sigma \in \mathcal{F}_{\text{int}}, \forall 1 \le i \le d. \qquad (6) $$

A network geometry is admissible if and only if it satisfies (2)–(4) and the family of measures is admissible. As soon as $\mathcal{N}$ is an admissible network, existence of an admissible network geometry was established in [12]. In [14] it is explained how to efficiently compute a network geometry by solving a symmetric linear system.

## 2. A NETWORK ELEMENT METHOD FOR ADVECTION-DIFFUSION-REACTION

In this section, after presenting our model problem we introduce the modified set of degrees of freedom for the modified NEM method that will allow to recover the conservation properties. We then propose modified versions of the original NEM reconstruction operators using those new degrees of freedom, and introduce a NEM reconstruction of the advective derivative to handle the advective terms. Then we propose a NEM scheme for advection-diffusion-reaction that takes as usual the form of a discrete variational problem mimicking the weak formulation of the model problem.

We consider the same diffusion-reaction model problem than in [13] complemented by a standard advection term:

$$ -\mathrm{div}(\Lambda \nabla u) + \mathrm{div}(\alpha u) + \beta u = f \quad \text{in } \Omega, \qquad (7) $$

with homogeneous Dirichlet boundary conditions $u = 0$ on $\partial\Omega$. The problem admits the following weak formulation, obtained from the most natural one through an integration by parts:

$$ \int_{\Omega} \Lambda \nabla u \nabla v + \frac{1}{2} \int_{\Omega} (\boldsymbol{\alpha} \cdot \nabla u) v - \frac{1}{2} \int_{\Omega} u (\boldsymbol{\alpha} \cdot \nabla v) + \int_{\Omega} \left( \frac{1}{2} \mathrm{div} \boldsymbol{\alpha} + \beta \right) uv = \int_{\Omega} fv \quad \forall v \text{ in } H_{0}^{1}(\Omega) $$
$$ \Leftrightarrow a(u, v) = a_{\Lambda}(u, v) + c_{\alpha}(u, v) + r_{\beta}(u, v) = l(v) \quad \forall v \text{ in } H_0^1(\Omega), \quad (8) $$

with obvious notations for the diffusion, advection and reaction bilinear forms $a_\Lambda$, $c_\alpha$ and $r_\beta$. We of course assume that $f \in L^2(\Omega)$, that $\beta$ is a measurable function from $\Omega$ to $\mathbb{R}$ with $0 \le \beta_* \le \beta(\mathbf{x}) \le \beta^*$ for almost every (a.e.) $\mathbf{x} \in \Omega$ where $\beta_*$, $\beta^*$ are two positive real numbers, that $\alpha \in H(\text{div}, \Omega)$ with $\text{div} \, \alpha \ge 0$ for a.e. $\mathbf{x} \in \Omega$, and that $\Lambda$ is a measurable function from $\Omega$ to $M_d(\mathbb{R})$ the set of $d \times d$ matrices and is such that for a.e. $\mathbf{x} \in \Omega$, $\Lambda(\mathbf{x})$ is symmetric positive definite and there exists two strictly positive real numbers $\lambda_*$, $\lambda^*$ such that for a.e. $\mathbf{x} \in \Omega$ and for every $\xi \in \mathbb{R}^d$ $\lambda_*|\xi|^2 \le \Lambda(\mathbf{x})\xi \cdot \xi \le \lambda^*|\xi|^2$.

Let $\mathcal{N}$ be an admissible discretization network endowed with an admissible network geometry $\mathcal{G}$. We associate to $\mathcal{N}$ the following space of degrees of freedom:

$$ X_{\mathcal{N}} = \{((u_K)_{K \in \mathcal{T}}, (u_\sigma)_{\sigma \in \mathcal{F}}) \mid u_K \in \mathbb{R}, \forall K \in \mathcal{T}, u_\sigma \in \mathbb{R}, \forall \sigma \in \mathcal{F}\}, $$

and denote $U = ((u_K)_{K \in \mathcal{T}}, (u_\sigma)_{\sigma \in \mathcal{F}})$. The difference with the original network element method of [12, 13] is the addition of the cell degrees of freedom $(u_K)_{K \in \mathcal{T}}$. Those degrees of freedom that will in some sense stand for the local mean value of the function will play a key role in deriving a conservative method. To account for homogeneous Dirichlet boundary conditions, we of course consider:

$$ X_{\mathcal{N},0} = \{U \in X_{\mathcal{N}} \mid u_\sigma = 0 \quad \text{for all } \sigma \in \mathcal{F}_{\text{ext}}\}, $$

and we define the local set of degrees of freedom associated to a cell as usual by setting:

$$ X_{\mathcal{N},K} = \{(u_K, (u_\sigma)_{\sigma \in \mathcal{F}_K}) \mid u_K \in \mathbb{R}, u_\sigma \in \mathbb{R}, \forall \sigma \in \mathcal{F}_K\} \quad \text{and} \quad X_{\mathcal{F}_K} = \{(u_\sigma)_{\sigma \in \mathcal{F}_K} \mid u_\sigma \in \mathbb{R}, \forall \sigma \in \mathcal{F}_K\}. $$

Of course, for any $U \in X_{\mathcal{N}}$, we denote $U_K = (u_K, (u_\sigma)_{\sigma \in \mathcal{F}_K})$ while for any $W \in X_{\mathcal{F}_K}^d$ we denote $W_{\mathcal{F}_K} = (w_\sigma)_{\sigma \in \mathcal{F}_K}$. The idea behind our modified NEM is to simply redefine the reconstruction operators of the original NEM and in particular the mean value operator $M_K$ taking advantage of the additional degree of freedom $u_K$ in each cell. Indeed, in each cell $K \in \mathcal{T}$ the mean value operator $M_K$ is now simply defined by setting

$$ \left| \begin{array}{l} \mathcal{M}_K : X_{\mathcal{N},K} \to \mathbb{P}_0(\mathbb{R}^d) \\ \mathbf{U}_K \to \mathcal{M}_K(\mathbf{U}_K) = u_K, \end{array} \right. \quad (9) $$

instead of using a barycentric interpolation from cell degrees of freedom as was done in [11–13]. The gradient operator is defined by:

$$ \left| \begin{array}{l} \nabla_K : X_{\mathcal{N},K} \to \mathbb{P}_0(\mathbb{R}^d)^d \\ \mathbf{U}_K \to \nabla_K(\mathbf{U}_K) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} (u_\sigma - \mathcal{M}_K(\mathbf{U}_K))\eta_{K,\sigma}. \end{array} \right. \quad (10) $$

Finally, the NEM local reconstruction operator $\Pi_K$ is defined for any cell $K \in \mathcal{T}$ as usual by:

$$ \left| \begin{array}{l} \Pi_K : X_{\mathcal{N},K} \to \mathbb{P}_1(\mathbb{R}^d) \\ \mathbf{U}_K \to \Pi_K(\mathbf{U}_K) = \mathcal{M}_K(\mathbf{U}_K) + \nabla_K(\mathbf{U}_K) \cdot (\mathbf{x} - \mathbf{x}_K), \end{array} \right. \quad (11) $$

implying that $\nabla \Pi_K(\mathbf{U}_K) = \nabla_K(\mathbf{U}_K)$. Given approximations $\alpha_{\mathcal{F}_K} = (\alpha_\sigma)_{\sigma \in \mathcal{F}_K}$ of $\alpha$ at the interfaces $\mathbf{x}_\sigma$ of $\mathcal{F}_K$, we define the advective derivative $\mathbf{G}_K$ by setting:

$$ \left| \begin{array}{l} \mathbf{G}_K : X_{\mathcal{F}_K}^d \times X_{\mathcal{N},K} \to \mathbb{P}_0(\mathbb{R}^d)^d \\ (\alpha_{\mathcal{F}_K}, \mathbf{U}_K) \to \mathbf{G}_K(\alpha_{\mathcal{F}_K}, \mathbf{U}_K) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \alpha_\sigma \cdot \eta_{K,\sigma} (u_\sigma - \mathcal{M}_K(\mathbf{U}_K)), \end{array} \right. \quad (12) $$
as well as the discrete divergence:

$$ \left| \begin{array}{l} \mathcal{D} \mathcal{I} \mathcal{V}_K : X_{\mathcal{F}_K}^d \to \mathbb{P}_0(\mathbb{R}^d)^d \\ \alpha_{\mathcal{F}_K} \to \mathcal{D} \mathcal{I} \mathcal{V}_K(\alpha_{\mathcal{F}_K}) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \alpha_\sigma \cdot \eta_{K,\sigma}. \end{array} \right. \qquad (13) $$

Of course for $\alpha$ regular enough we will use $\alpha_\sigma = \alpha(x_\sigma)$ in practice. The introduction of the operator $M_K$ that is simply defined as equal to the degree of freedom $u_K$ has the advantage of establishing a formal equivalence between the modified NEM and the original one. In practice, it means that many proofs that were conducted using $M_K$ for the original version can be directly transposed to the modified version.

**Remark 2.1.** In the non conservative case [12], the local gradient operator was defined by setting:

$$ \nabla_K (U_K) = \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} u_\sigma \eta_{K,\sigma}. \qquad (14) $$

Using the original mean value operator $M_K$ defined through barycentric interpolation, it is interesting to notice that one could use (10) instead of (14) also in the original setting, and all the results of [12] would still hold including the approximation order. In fact, in the case of an exact geometry, as

$$ \sum_{\sigma \in \mathcal{F}_K} \eta_{K,\sigma} = 0, $$

the two definitions coincide. It is only when the geometry is only approximate that a geometric consistency error separates the two definitions.

The NEM approximation of (8) consists in finding a solution $U \in X_{N,0}$ of

$$ a_{\Lambda,h}(U, V) + c_{\alpha,h}(U, V) + r_{\beta,h}(U, V) = l_h(V) \quad \text{for all } V \in X_{N,0}. \qquad (15) $$

The approximation of the diffusion-reaction $a_{\Lambda,h} + r_{\beta,h}$ will be formally identical to the one of the original NEM of [13] simply replacing $M_K$ and $\nabla \Pi_K$ by their new definitions involving the degrees of freedom $(u_K)_{K \in \mathcal{T}}$, where we recall that (see [13]):

$$ a_{\Lambda,h}(U, V) = \sum_{K \in \mathcal{T}} a_{\Lambda,h}^K(U_K, V_K) \quad \text{and} \quad r_{\beta,h}(U, V) = \sum_{K \in \mathcal{T}} r_h^K(U_K, V_K), $$

with

$$ a_{\Lambda,h}^K(U_K, V_K) = m_K \Lambda_K \nabla \Pi_K(U_K) \cdot \nabla \Pi_K(V_K) + s^K(U_K - \mathcal{D}_K(\Pi_K(U_K)), V_K - \mathcal{D}_K(\Pi_K(V_K))), \quad (16) $$

and

$$ r_{\beta,h}^K(U_K, V_K) = m_K \beta_K M_K(U_K) M_K(V_K) + s_0^K(U_K - \mathcal{D}_K(M_K(U_K)), V_K - \mathcal{D}_K(M_K(V_K))), \quad (17) $$

where $\Lambda_K$ is an approximation of $\Lambda$ around $x_K$ and $\beta_K$ is an approximation of $\beta$ around $x_K$, and where the stabilizations $s^K$ and $s_0^K$ are positive symmetric bilinear forms on $X_{N,K} \times X_{N,K}$, such that

$$ s^K(U_K, V_K) = m_K h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} S_{K,\sigma,\sigma'} u_\sigma v_{\sigma'} \quad \text{and} \quad s_0^K(U_K, V_K) = m_K \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} S_{K,\sigma,\sigma'}^0 u_\sigma v_{\sigma'}, \quad (18) $$
where $S_K = (S_{K,\sigma,\sigma'})_{\sigma,\sigma' \in \mathcal{F}_K}$ (resp. $S_K^0 = (S_{K,\sigma,\sigma'}^0)_{\sigma,\sigma' \in \mathcal{F}_K}$) can be any symmetric positive definite matrix independent on the geometry $\mathcal{G}$ associated to the network. For the right-hand side, assuming that $f_K$ is again an approximation of $f$ at $x_K$ we set

$$l_h(\mathbf{V}) = \sum_{K \in \mathcal{T}} m_K f_K \mathcal{M}_K(\mathbf{V}_K).$$

This implies that all the results valid for the original NEM and that did not involve the explicit definitions of $\mathcal{M}_K$ and $\nabla \Pi_K$ will remain valid in the present case. In particular, the stability results for the diffusion-reaction bilinear form directly carry over to the modified version.

The novelty in the NEM approximation is the advection bilinear form $c_{\alpha,h}$ defined by:

$$c_{\alpha,h}(\mathbf{U}, \mathbf{V}) = \sum_{K \in \mathcal{T}} c_{\alpha,h}^{K}(\mathbf{U}_K, \mathbf{V}_K),$$

with $c_{\alpha,h}^K$ mimicking $c_\alpha$ using the discrete advective derivative $\mathcal{G}_K$ and the discrete divergence $\mathcal{D}\mathcal{I}\mathcal{V}_K$:

$$c_{\alpha,h}^{K}(\mathbf{U}_K, \mathbf{V}_K) = \frac{1}{2} m_K G(\alpha_{\mathcal{F}_K}, \mathbf{U}_K) \mathcal{M}_K(\mathbf{V}) - \frac{1}{2} m_K G(\alpha_{\mathcal{F}_K}, \mathbf{V}_K) \mathcal{M}_K(\mathbf{U}) \\ + \frac{1}{2} m_K (\mathcal{D}\mathcal{I}\mathcal{V}_K(\alpha_{\mathcal{F}_K}))^+ \mathcal{M}_K(\mathbf{U}_K) \mathcal{M}_K(\mathbf{V}_K) + \frac{1}{2} \sum_{\sigma \in \mathcal{F}_K} |\alpha_\sigma \cdot \eta_{K,\sigma}| (u_\sigma - \mathcal{M}_K(\mathbf{U}))(v_\sigma - \mathcal{M}_K(\mathbf{V})),$$

where $x^+ = \max(0, x)$. The last term is a stabilization bilinear form whose role will appear more obvious when rewriting the scheme under flux form. Notice that thanks to the formal analogy with the original NEM reconstruction operators, using the original definitions of NEM operators the above bilinear form $c_{\alpha,h}$ also explains how to approximate advection terms for the original NEM using the definition of $\mathcal{M}_K$ of [12, 13].

### 3. SOME REMARKS ON STABILITY AND CONVERGENCE

The objective of this section is to briefly explain why the above NEM method is stable and which error estimates can be obtained, relying on the formal analogy with the original NEM method for diffusion reaction. Recalling the definition of the NEM local semi-norm:

$$\|\mathbf{U}\|_{1,K}^2 = \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} |u_\sigma - \mathcal{M}_K(\mathbf{U}_K)|^2,$$

thanks to the analogy with non-conservative version we know from [12] that there exists $\gamma > 0$ depending $\Lambda, \beta$ and the quality parameters of the network but not on $h$ such that:

$$a_{\Lambda,h}^{K}(\mathbf{U}_K, \mathbf{U}_K) \geq \gamma \|\mathbf{U}\|_{1,K}^2.$$

Obviously, we have $r_\beta(\mathbf{U}_K, \mathbf{U}_K) \geq 0$, and as $c_{\alpha,h}^K(\mathbf{U}_K, \mathbf{U}_K)$ reduces to:

$$c_{\alpha,h}^{K}(\mathbf{U}_K, \mathbf{U}_K) = \frac{1}{2} m_K (\mathcal{D}\mathcal{I}\mathcal{V}_K(\alpha_{\mathcal{F}_K}))^+ |\mathcal{M}_K(\mathbf{U}_K)|^2 + \frac{1}{2} \sum_{\sigma \in \mathcal{F}_K} |\alpha_\sigma \cdot \eta_{K,\sigma}| |u_\sigma - u_K|^2 \geq 0,$$

it is clear that $a_h^K(\mathbf{U}_K, \mathbf{U}_K) = a_{\Lambda,h}^K(\mathbf{U}_K, \mathbf{U}_K) + c_{\alpha,h}^K(\mathbf{U}_K, \mathbf{U}_K) + r_{\beta,h}^K(\mathbf{U}_K, \mathbf{U}_K) \geq \gamma \|\mathbf{U}\|_{1,K}^2$, ensuring that the method is unconditionally stable.

We know consider the approximation properties of the new NEM method and focus on the advection bilinear form. For any $\varphi \in C^0(\mathbb{R}^d)$ and any $\Phi \in C^0(\mathbb{R}^d)^d$ (and more generally for any function for which it makes sense), we denote

$$\mathcal{D}_K(\varphi) = (\varphi(\mathbf{x}_K), (\varphi(\mathbf{x}_\sigma))_{\sigma \in \mathcal{F}_K}) \quad \text{and} \quad \mathcal{D}_{\mathcal{F}_K}(\Phi) = (\Phi(\mathbf{x}_\sigma))_{\sigma \in \mathcal{F}_K},$$
the local sets of degrees of freedom associated with $\varphi$ and $\Phi$. We know from [12] that thanks to (2)-(3) the discrete divergence operator is consistent, i.e. for instance if $\alpha \in C^2(\Omega)^d$ there exists $C_\alpha > 0$ such that for all $\mathbf{x} \in B(\mathbf{x}_K, h_K)$:

$$|\mathcal{D}\mathcal{I}\mathcal{V}_K(\mathcal{D}_{\mathcal{F}_K}(\alpha)) - \text{div } \alpha(\mathbf{x})| \le C_\alpha h_K.$$

For div $\alpha(\mathbf{x}) \ge 0$, as $|\mathcal{D}\mathcal{I}\mathcal{V}_K(\mathcal{D}_{\mathcal{F}_K}(\alpha))^-| \le |\mathcal{D}\mathcal{I}\mathcal{V}_K(\mathcal{D}_{\mathcal{F}_K}(\alpha)) - \text{div } \alpha(\mathbf{x})|$ it is clear that:

$$|\mathcal{D}\mathcal{I}\mathcal{V}_K(\mathcal{D}_{\mathcal{F}_K}(\alpha))^+ - \text{div } \alpha(\mathbf{x})| = |\mathcal{D}\mathcal{I}\mathcal{V}_K(\mathcal{D}_{\mathcal{F}_K}(\alpha)) + \mathcal{D}\mathcal{I}\mathcal{V}_K(\mathcal{D}_{\mathcal{F}_K}(\alpha))^- - \text{div } \alpha(\mathbf{x})| \le 2C_\alpha h_K.$$

Next, assume that both $\alpha$ and $u$ are first order polynomials. Then:

$$
\begin{align*}
\mathcal{G}_K(\mathcal{D}_{\mathcal{F}_K}(\alpha), \mathcal{D}_K(u)) &= \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} \alpha(\mathbf{x}_\sigma) \cdot \eta_{K,\sigma}(u(\mathbf{x}_\sigma) - u(\mathbf{x}_K)) \\
&= \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} (\alpha(\mathbf{x}_K) + \nabla\alpha(\mathbf{x}_K)(\mathbf{x}_\sigma - \mathbf{x}_K)) \cdot \eta_{K,\sigma}(\nabla u(\mathbf{x}_K) \cdot (\mathbf{x}_\sigma - \mathbf{x}_K)) \\
&= \frac{1}{m_K} \alpha(\mathbf{x}_K) \cdot \left( \sum_{\sigma \in \mathcal{F}_K} \eta_{K,\sigma}(\mathbf{x}_\sigma - \mathbf{x}_K)^T \right) \nabla u(\mathbf{x}_K) \\
&\quad + \frac{1}{m_K} \sum_{\sigma \in \mathcal{F}_K} (\nabla\alpha(\mathbf{x}_K)(\mathbf{x}_\sigma - \mathbf{x}_K)) \cdot \eta_{K,\sigma}(\nabla u(\mathbf{x}_K) \cdot (\mathbf{x}_\sigma - \mathbf{x}_K)).
\end{align*}
$$

From (3) as the second term is clearly bounded by $h_K^2$ times a constant depending only on $u$, $\alpha$ and the quality parameter $\theta_\Pi = \sup_{K \in \mathcal{T}} \sup_{\sigma \in \mathcal{F}_K} m_K^{-1} h_K |\eta_{K,\sigma}|$ we clearly get $\mathcal{G}_K(\mathcal{D}_{\mathcal{F}_K}(\alpha), \mathcal{D}_K(u)) \approx \alpha(\mathbf{x}_K) \cdot \nabla u(\mathbf{x}_K)$ as expected. Building on this observation, using Taylor expansion and proceeding as in [12], it is straightforward and classical to establish the following strong consistency result for the advective derivative:

**Proposition 3.1.** Let $(\mathcal{N}, \mathcal{G})$ be an admissible network and associated admissible geometry. Then for any $\alpha \in C^2(\bar{\Omega})^d$ and any $u \in C^2(\bar{\Omega})$ there exists $C > 0$ depending only on $\theta_\Lambda$, $\theta_\Pi$, $\alpha$ and $u$ such that for any $\mathbf{x} \in B(\mathbf{x}_K, h_K)$

$$|\mathcal{G}_K(\mathcal{D}_{\mathcal{F}_K}(\alpha), \mathcal{D}_K(u)) - \alpha(\mathbf{x}) \cdot \nabla u(\mathbf{x})| \le C h_K,$$

with obvious extension to cases with Sobolev regularity only along the lines of [11]. To develop a full convergence theory for the conservative network element method, it remains to follow the lines of the non conservative case [11, 13] for diffusion reaction, using the strong consistency of the additional advective terms emphasized above. In fact, once again by virtue of the formal equivalence between the two schemes, a comprehensive convergence study would only differ from the non conservative case through slight technical differences in the proof of the estimates, due to the differences in the definitions of the operators $\mathcal{M}_K$ and $\Pi_K$. In particular, it should be obvious from [12] that the newly defined operator $\Pi_K$ shares the strong and weak consistency properties of the original one, while the new $\mathcal{M}_K$ is obviously strongly consistent. Elaborating along the lines of [11, 13], one ultimately obtains convergence to minimal regularity solutions and the expected error estimates for the network element interpolants defined in [11]. To give some highlights on the result, we mention in particular that we can obtain the following local estimates as a corollary of a general error estimate:

**Proposition 3.2.** Let $(\mathcal{N}, \mathcal{G})$ be an admissible discretization network and an associated admissible network geometry, and let $\mathcal{U}$ be the solution of the associated problem (15). Assume that there exists $C_f > 0$, $C_\Lambda > 0$, $C_\beta > 0$ and $C_\alpha > 0$ such that:

$$\left( \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} |f_K - f|^2 \right)^{1/2} \le C_f h, \quad \sup_{K \in \mathcal{T}} \|\beta - \beta_K\|_{L^\infty(B_K)} \le C_\beta h \quad \text{and} \quad \sup_{K \in \mathcal{T}} \|\Lambda - \Lambda_K\|_{L^\infty(B_K)} \le C_\Lambda h,$$
as well as

$$ \sup_{\sigma \in \mathcal{F}} \sup_{K \in \mathcal{T}_{\sigma}} \|\boldsymbol{\alpha}_{\sigma} - \boldsymbol{\alpha}\|_{L^{\infty}(B_K)} \le C_{\alpha} h \quad \text{and} \quad \sup_{K \in \mathcal{T}} \|\mathcal{D} \mathcal{I} \mathcal{V}_K(\boldsymbol{\alpha}_{\mathcal{F}_K}) - \text{div} \boldsymbol{\alpha}\|_{L^{\infty}(B_K)} \le C_{\alpha} (h + h^p), $$

and assume that the solution $u$ of (8), $\alpha$ and $\Lambda$ are regular enough. Then, for any $K \in \mathcal{T}$

$$ \|u - \mathcal{M}_K(\mathbf{U})\|_{L^2(B_K \cap \Omega)} \le C(h + h^p) \quad \text{and} \quad \|\nabla u - \nabla_K(\mathbf{U})\|_{L^2(B_K \cap \Omega)^d} \le C(h + h^p), \quad (19) $$

where the constant $C > 0$ depends on $\Omega$, $u$, $\alpha$, $\Lambda$, $\beta$ and the quality parameters of the network but not on $h$.

As it would otherwise involve a long exposition of the network element convergence framework, to get a precise statement on the needed regularity we refer the reader to [13] and to Appendix B where we sketch the proof of this result. Notice that it heavily relies on the coercivity of the diffusion operator. Consequently, the constants in the above estimates strongly depends on the diffusion tensor $\Lambda$. For vanishing diffusion, it is well known that a loss of convergence order should be expected if the advection is not zero (see [16] for instance). Similarly to what can be done for mesh-based approaches, we expect that convergence estimates robust to vanishing diffusion relying on the discrete version of the advective derivative $\boldsymbol{\alpha} \cdot \nabla$ could be derived for the NEM. As our goal here is essentially to present the modified NEM scheme and as it would otherwise forces us to a lengthy derivation of convergence estimates in an advection adapted discrete norm, we will only explore this limit behavior through numerical experiments.

## 4. A FLUX FORMULATION FOR THE NEW NETWORK ELEMENT METHOD

The aim of this section is to explain how to derive a flux formulation for the modified NEM, that will now justify the “conservative” claim for the new method. This flux formulation will both be much more familiar for finite volume practitioners and provide an easy way for implementing the conservative NEM. In particular, we will see that the advective term almost reduces to the usual upwind flux term, thanks to the presence of the stabilization in the bilinear form $c_{\alpha,h}^K$.

We start by the diffusive term. For any $(\mathbf{U}_K, \mathbf{V}_K) \in X_{\mathcal{N},K} \times X_{\mathcal{N},K}$ we have:

$$ a_{\Lambda,h}^{K}(\mathbf{U}_K, \mathbf{V}_K) = m_K \Lambda_K \nabla \Pi_K(\mathbf{U}) \nabla \Pi_K(\mathbf{V}) + m_K h_K^{-2} \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} S_{K,\sigma,\sigma'}(u_\sigma - \Pi_K(\mathbf{U})(\mathbf{x}_\sigma))(v_{\sigma'} - \Pi_K(\mathbf{V})(\mathbf{x}_\sigma)). $$

Notice that by definition of $\Pi_K$, we have:

$$ u_{\sigma} - \Pi_{K}(\mathbf{U}_{K})(\mathbf{x}_{\sigma}) = u_{\sigma} - u_{K} - \frac{1}{m_{K}} \sum_{\sigma' \in \mathcal{F}_{K}} (u_{\sigma'} - u_{K}) \eta_{K, \sigma'} \cdot (\mathbf{x}_{\sigma} - \mathbf{x}_{K}). $$

Consequently it is clear that expanding everything, we will get an expression that takes the form:

$$ a_{\Lambda,h}^{K}(\mathbf{U}_K, \mathbf{V}_K) = \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} \mathcal{A}_K^{\sigma, \sigma'} (u_{\sigma} - u_{K})(v_{\sigma'} - v_{K}), $$

where only remains terms in $u_\sigma - u_K$ and $v_{\sigma'} - v_K$. This is precisely what is needed to obtain a flux formulation. More precisely, a straightforward computation leads to:

$$ \left\{
\begin{aligned}
\mathcal{A}_K^{\sigma, \sigma'} &= \frac{1}{m_K} \Lambda_K \eta_{K, \sigma} \cdot \eta_{K, \sigma'} + m_K h_K^{-2} \mathcal{S}_K^{\sigma, \sigma'}, \\
\mathcal{S}_K^{\sigma, \sigma'} &= y_{K, \sigma'}^T S_K y_{K, \sigma}, \\
y_{K, \sigma} &= (y_{K, \sigma}^{\prime\prime})_{\sigma'' \in \mathcal{F}_K}, \quad y_{K, \sigma}^{\prime\prime} = \delta_{\sigma, \sigma''} - \frac{1}{m_K} \eta_{K, \sigma} \cdot (\mathbf{x}_{\sigma''} - \mathbf{x}_K),
\end{aligned}
\right.
\qquad (20) $$
where $\delta_{\sigma,\sigma''}$ is the Kronecker symbol. The diffusion fluxes are then defined by setting

$$F_{K,\sigma}^{D}(\mathbf{U}_K) = \sum_{\sigma' \in \mathcal{F}_K} \mathcal{A}_K^{\sigma',\sigma}(u_K - u_{\sigma'}), \qquad (21)$$

in order to finally obtain:

$$a_{A,h}^{K}(\mathbf{U}_K, \mathbf{V}_K) = \sum_{\sigma \in \mathcal{F}_K} F_{K,\sigma}^{D}(\mathbf{U})(v_K - v_{\sigma}). \qquad (22)$$

In the same way, we get for any $(\mathbf{U}_K, \mathbf{V}_K) \in X_{N,K} \times X_{N,K}$:

$$r_{\beta,h}^{K}(\mathbf{U}_K, \mathbf{V}_K) = m_K \beta_K u_K v_K + m_K \sum_{\sigma \in \mathcal{F}_K} \sum_{\sigma' \in \mathcal{F}_K} S_{K,\sigma,\sigma'}^{0}(u_{\sigma} - u_{K})(v_{\sigma'} - v_K).$$

After a straightforward computation, this leads to:

$$r_{\beta,h}^{K}(\mathbf{U}_K, \mathbf{V}_K) = m_K \beta_K u_K v_K + \sum_{\sigma \in \mathcal{F}_K} F_{K,\sigma}^{R}(\mathbf{U})(v_K - v_{\sigma}),$$

where the "reaction fluxes" are defined by:

$$F_{K,\sigma}^{R}(\mathbf{U}_K) = m_K \sum_{\sigma' \in \mathcal{F}_K} S_{K,\sigma',\sigma}^{0}(u_K - u_{\sigma'}). \qquad (23)$$

Now we turn to the more interesting and involved advective terms. Expanding the various terms, we get for any $(\mathbf{U}_K, \mathbf{V}_K) \in X_{N,K} \times X_{N,K}$:

$$c_{\alpha,h}^{K}(\mathbf{U}_K, \mathbf{V}_K) = \frac{1}{2} \sum_{\sigma \in \mathcal{F}_K} \boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma}(u_{\sigma} - u_K)v_K - \frac{1}{2} \sum_{\sigma \in \mathcal{F}_K} \boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma}(v_{\sigma} - v_K)u_K + \frac{1}{2} \left( \sum_{\sigma \in \mathcal{F}_K} \boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} \right) u_K v_K \\ + \frac{1}{2} m_K \left( (\mathcal{D}\mathcal{I}\mathcal{V}(\boldsymbol{\alpha}_{\mathcal{F}_K}))^{+} - \mathcal{D}\mathcal{I}\mathcal{V}(\boldsymbol{\alpha}_{\mathcal{F}_K}) \right) u_K v_K + \frac{1}{2} \sum_{\sigma \in \mathcal{F}_K} |\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma}| (u_{\sigma} - u_K)(v_{\sigma} - v_K),$$

where we have used $(\mathcal{D}\mathcal{I}\mathcal{V}(\boldsymbol{\alpha}_{\mathcal{F}_K}))^{+} = \mathcal{D}\mathcal{I}\mathcal{V}(\boldsymbol{\alpha}_{\mathcal{F}_K}) + (\mathcal{D}\mathcal{I}\mathcal{V}(\boldsymbol{\alpha}_{\mathcal{F}_K}))^{+} - \mathcal{D}\mathcal{I}\mathcal{V}(\boldsymbol{\alpha}_{\mathcal{F}_K})$. We now regroup the terms in front of $v_K$ and those in front on $v_{\sigma}$, leading to:

$$c_{\alpha,h}^{K}(\mathbf{U}_K, \mathbf{V}_K) = \frac{1}{2} m_K ((\mathcal{D}\mathcal{I}\mathcal{V}(\boldsymbol{\alpha}_{\mathcal{F}_K}))^{+} - \mathcal{D}\mathcal{I}\mathcal{V}(\boldsymbol{\alpha}_{\mathcal{F}_K})) u_K v_K \\ + \frac{1}{2} \sum_{\sigma \in \mathcal{F}_K} ((\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} + |\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma}|) u_K + (\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} - |\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma}|) u_{\sigma}) v_K \\ - \frac{1}{2} \sum_{\sigma \in \mathcal{F}_K} ((\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} - |\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma}|) u_{\sigma} + (\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} + |\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma}|) u_K) v_{\sigma} + \frac{1}{2} \sum_{\sigma \in \mathcal{F}_K} \boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} u_{\sigma} v_{\sigma}.$$

The advective flux are defined by setting:

$$F_{K,\sigma}^{A}(\mathbf{U}_K) = \frac{1}{2} (\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} + |\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma}|) u_K + (\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} - |\boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma}|) u_{\sigma},$$

which rewrites under the more familiar upwind form:

$$F_{K,\sigma}^{A}(\mathbf{U}_K) = \boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} u_{\sigma}^{\uparrow} \quad \text{where} \quad u_{\sigma}^{\uparrow} = \begin{cases} u_K & \text{if } \boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} \ge 0, \\ u_{\sigma} & \text{if } \boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} < 0. \end{cases}$$
Notice that we obtain the familiar upwind scheme only thanks to the presence of the stabilization term in the bilinear form $c_{\alpha,h}^K$, which is precisely the reason why we have introduced this term. Let us finally denote:

$$ \delta_K(\boldsymbol{\alpha}) = \frac{1}{2}(\mathcal{D}\mathcal{I}\mathcal{V}(\boldsymbol{\alpha}_{\mathcal{F}_K}))^+ - \mathcal{D}\mathcal{I}\mathcal{V}(\boldsymbol{\alpha}_{\mathcal{F}_K})) = \frac{1}{2m_K} \left( \left( \sum_{\sigma \in \mathcal{F}_K} \boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} \right)^+ - \sum_{\sigma \in \mathcal{F}_K} \boldsymbol{\alpha}_{\sigma} \cdot \boldsymbol{\eta}_{K,\sigma} \right). $$

Gathering those results, the discrete variational formulation (15) becomes:

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} (F_{K,\sigma}^D(U_K) + F_{K,\sigma}^R(U_K) + F_{K,\sigma}^A(U_K))(v_K - v_\sigma) + \sum_{K \in \mathcal{T}} m_K(\beta_K + \delta_K(\boldsymbol{\alpha}))u_K v_K \\ = \sum_{K \in \mathcal{T}} m_K f_K v_K - \frac{1}{2} \sum_{\sigma \in \mathcal{F}} \boldsymbol{\alpha}_{\sigma} \cdot \left( \sum_{K \in \mathcal{T}_{\sigma}} \boldsymbol{\eta}_{K,\sigma} \right) u_{\sigma} v_{\sigma}. $$

Taking one degree of freedom of $\mathcal{V}$ equal to one and all others equal to 0 and doing so for every dof, we obtain:

$$ \left| \begin{array}{l} m_K(\beta_K + \delta_K(\boldsymbol{\alpha}))u_K + \sum_{\sigma \in \mathcal{F}_K} F_{K,\sigma}^D(U_K) + F_{K,\sigma}^R(U_K) + F_{K,\sigma}^A(U_K) = m_K f_K \quad \forall K \in \mathcal{T}, \\ \sum_{K \in \mathcal{T}_{\sigma}} F_{\sigma,K}^D(U_K) + F_{\sigma,K}^R(U_K) + F_{\sigma,K}^A(U_K) = -\frac{1}{2} \boldsymbol{\alpha}_{\sigma} \cdot \left( \sum_{K \in \mathcal{T}_{\sigma}} \boldsymbol{\eta}_{K,\sigma} \right) u_{\sigma} \quad \forall \sigma \in \mathcal{F}_{\text{int}}, \\ u_{\sigma} = 0 \quad \forall \sigma \in \mathcal{F}_{\text{ext}}, \end{array} \right. \quad (24) $$

where we have of course denoted $F_{\sigma,K}^A(U_K) = -F_{K,\sigma}^A(U_K)$, $F_{\sigma,K}^D(U_K) = -F_{K,\sigma}^D(U_K)$ and $F_{\sigma,K}^R(U_K) = -F_{K,\sigma}^R(U_K)$. In other words our method shares the same analytic formulation than classical finite volume methods, with the cells exchanging with the interfaces through fluxes and reciprocally. We also recover the local balance of diffusion-reaction fluxes ensured by $F_{\sigma,K}^D = -F_{K,\sigma}^D$ and $F_{\sigma,K}^R = -F_{K,\sigma}^R$, which immediately implies global conservation for the pure diffusion-reaction case. Indeed assume that $\boldsymbol{\alpha} = 0$, then summing over $K \in \mathcal{T}$ and adding the sum over $\sigma \in \mathcal{F}_{\text{int}}$ of the flux continuity equations we get:

$$ \sum_{K \in \mathcal{T}} m_K \beta_K u_K + \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} F_{K,\sigma}^D(U_K) + F_{K,\sigma}^R(U_K) = \sum_{K \in \mathcal{T}} m_K f_K, $$

which justifies the claim of a conservative network element method, at least for pure diffusion-reaction problems. In the presence of advection, proceeding in the same way leads to:

$$ \sum_{K \in \mathcal{T}} m_K(\beta_K + \delta_K(\boldsymbol{\alpha}))u_K + \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_{\text{ext}}} F_{K,\sigma}^D(U_K) + F_{K,\sigma}^R(U_K) + F_{K,\sigma}^A(U_K) \\ = \sum_{K \in \mathcal{T}} m_K f_K - \frac{1}{2} \sum_{\sigma \in \mathcal{F}} \boldsymbol{\alpha}_{\sigma} \cdot \left( \sum_{K \in \mathcal{T}_{\sigma}} \boldsymbol{\eta}_{K,\sigma} \right) u_{\sigma}. $$

For admissible only geometry, the global conservation error induced by the approximate compatibility is controlled by $\theta_A h^p$, and will thus rapidly disappear with an improvement in the discrete geometry. For exactly compatible geometries, i.e. such that the compatibility errors $\epsilon_\sigma = 0$, we have $\sum_{K \in \mathcal{T}_\sigma} \boldsymbol{\eta}_{K,\sigma} = 0$ and thus we completely get rid of this term. Fortunately, from [14] we know that one can always compute an exactly compatible geometry (up to solver precision). We nevertheless found interesting to derive everything in the most general case and also to notice that for pure diffusion-reaction, even without exact compatibility global conservation is achieved.
Another contribution to the global conservation error comes from the additional numerical reaction coefficient $\delta_K(\alpha)$. This term is not linked to the geometry error $\epsilon(\mathcal{G})$ but rather comes from the fact that contrary to mesh-based methods, our discrete divergence operator applied to the degrees of freedom of $\alpha$ do not coincide with the $L^2$ projection on cell-wise constant vectors of the exact divergence. Thus to rigorously ensure the stability of the NEM, we had to incorporate this term to account for the potential negative values of $\mathcal{DIV}(\alpha_{\mathcal{F}_K})$. For non-vanishing diffusion and/or reaction coefficients, unless the mesh size is high in general one can simply use $\delta_K(\alpha) = 0$ and still obtain a stable method, which amounts to directly using $\mathcal{DIV}_K(\alpha_{\mathcal{F}_K})$ instead of $|\mathcal{DIV}_K(\alpha_{\mathcal{F}_K})|$ in the advection bilinear form.

**Remark 4.1.** It is of course possible to replace the term $|\mathcal{DIV}_K(\alpha_{\mathcal{F}_K})|$ by some approximation $(\text{div } \alpha)_K$ taking positive values in the advection bilinear form. This ultimately simply leads to modify the value $\delta_K(\alpha)$ used in the finite volume like formulation by $\delta_K(\alpha) = \frac{1}{2m_K}((\text{div } \alpha)_K - \mathcal{DIV}(\alpha_{\mathcal{F}_K}))$. Also notice that all the results of the present paper remain of course valid under the more general hypothesis $\frac{1}{2}\text{div } \alpha + \beta \ge 0$. In this case, the term $\delta_K(\alpha)$ should be replaced by

$$ \delta_K(\alpha, \beta) = \left( \left( \frac{1}{2} \mathcal{DIV}(\alpha_{\mathcal{F}_K}) + \beta_K \right)^+ - \left( \frac{1}{2} \mathcal{DIV}(\alpha_{\mathcal{F}_K}) + \beta_K \right) \right). $$

Within our notations, this implies mixing the advective bilinear form $c_{\alpha,h}^K$ with the reaction coefficient $\beta_K$ in principle handled by the reactive bilinear form $r_{\beta,h}^K$. Although this is more general, we felt that it would ease the presentation to be able to fully separate the discretization of advection and reaction, hence the more restrictive hypothesis $\text{div } \alpha \ge 0$.

Apart from exact compatibility of the geometry, the key to obtain a conservative NEM was to introduce the additional cell degrees of freedom $(u_K)_{K \in \mathcal{T}}$ that plays the role of Lagrange multipliers for the local flux balance property. Notice that in practice, the additional cell unknowns $(u_K)_{K \in \mathcal{T}}$ can be locally eliminated through static condensation (i.e. local inversion of the cell equations) thus reducing the system to one that has the same size as the original network element system. Recovering the conservation property thus only implies a minimal computational overhead. The general idea of the implementation of the method is expressed through Algorithm 1: it exactly matches the usual implementation of the Hybrid finite volume scheme (see [22]).

## 5. A CONSERVATIVE AND STABLE PROJECTION BETWEEN NETWORKS

Before turning to numerical experiments, let us explain how to perform efficiently a conservative projection between networks. This will indeed be a necessary step when considering Eulerian simulations on moving domains. Let $(\mathcal{N}, \mathcal{G})$ and $(\mathcal{N}, \hat{\mathcal{G}})$ be two pairs of admissible network and associated geometry. To project an element $\hat{U}$ of $X_{\mathcal{N}}$ defined on $(\mathcal{N}, \hat{\mathcal{G}})$ into an element $P_{\mathcal{N} \to \mathcal{N}}(\hat{U})$ of $X_{\mathcal{T}} = \{(u_K)_{K \in \mathcal{T}} | u_K \in \mathbb{R} | K \}$ defined on $(\mathcal{N}, \mathcal{G})$, we consider the following minimization problem:

$$ P_{\mathcal{N} \to \mathcal{N}}(\hat{U}) = \arg \min_{U \in C_{\mathcal{T}}(\hat{U})} J_{\mathcal{N} \to \mathcal{N}}(U, \hat{U}), \qquad (25) $$

where the cost function $J_{\mathcal{N} \to \mathcal{N}}$ is defined by:

$$ J_{\mathcal{N} \to \mathcal{N}}(U, \hat{U}) = \frac{1}{2} \sum_{K \in \mathcal{T}} m_K \left( u_K - \Pi_{\mathcal{N},K}(\hat{U})(x_K) \right)^2, \qquad (26) $$

and the constrained set $C_{\mathcal{T}}(\hat{U})$ correspond to

$$ C_{\mathcal{T}}(\hat{U}) = \left\{ \mathbf{V} \in X_{\mathcal{T}} \mid \sum_{K \in \mathcal{T}} m_K v_K = \sum_{\hat{K} \in \hat{\mathcal{T}}} m_{\hat{K}} \hat{u}_{\hat{K}} \quad \text{and} \quad \sum_{K \in \mathcal{T}} m_K v_K^2 = \sum_{\hat{K} \in \hat{\mathcal{T}}} m_{\hat{K}} \hat{u}_{\hat{K}}^2 \right\}. \qquad (27) $$
**Algorithm 1** Pseudo code for the conservative network element method

1: Initialization:

2: All matrix and vector coefficients are set to zero.

3: Assembly loop:

4: **for** K ∈ T **do**

5: Compute $\delta_K(\alpha)$

6: Add the cell terms $m_K(\beta_K + \delta_K(\alpha))u_K$ and second member $m_K f_K$ to the cell equations of (24).

7: **for** $\sigma, \sigma' \in \mathcal{F}_K$ **do**

8: Compute all the local terms $A_K^{\sigma, \sigma'}$.

9: **end for**

10: **for** $\sigma \in \mathcal{F}_K$ **do**

11: Form the local flux $F_{\sigma,K}^A(U_K)$, $F_{\sigma,K}^R(U_K)$ and $F_{\sigma,K}^D(U_K)$ and add their contributions to the cell equation on K and the face equation on $\sigma$ of (24).

12: **end for**

13: **end for**

14: **for** $\sigma \in \mathcal{F}_{\text{int}}$ **do**

15: Add the face term $-\frac{1}{2}\alpha_{\sigma} \cdot (\sum_{K \in \mathcal{T}_{\sigma}} \eta_{K,\sigma}) u_{\sigma}$ to the face equation on $\sigma$.

16: **end for**

17: **for** $\sigma \in \mathcal{F}_{\text{ext}}$ **do**

18: Set the face equation on $\sigma$ to be the Dirichlet boundary condition equation.

19: **end for**

20: Static condensation:

21: **for** K ∈ T **do**

22: Solve the cell equations of (24) to express $u_K$ from $(u_{\sigma})_{\sigma \in \mathcal{F}_K}$.

23: **end for**

24: Linear solver step: solve the reduced linear system formed by the interface equations of (24).

In other words, we seek an approximation of $\hat{U}$ in $X_T$ that ensures global conservativity and stability in $L^2$ discrete norm. The quality (or consistency) of the projection will be ensured through the operator $\Pi_{\hat{K},K}$. Many choices are possible for this operator, in the present paper we consider the “nearest neighbour” choice corresponding to using the second order reconstruction operator $\Pi_{\hat{K}}$ associated to the cell $\hat{K} \in \hat{T}$ that minimizes the distance $|x_{\hat{K}} - x_K|$, i.e.:

$$ \Pi_{\hat{K},K}(\hat{U}) = \Pi_{\hat{K}(K)}(\hat{U}) \quad \text{where } \hat{K}(K) \in \hat{T} \text{ is such that } |x_{\hat{K}(K)} - x_K| = \min_{\hat{K} \in \hat{T}} |x_{\hat{K}} - x_K|. \quad (28) $$

Another interesting choice could be to use a weighted sum of the second order approximations $\Pi_{\hat{K}}(\hat{U})(x_K)$ for all $\hat{K} \in \hat{T}$ such that $B_{\hat{K}} \cap B_K \neq 0$, with weights related to powers of $|x_{\hat{K}} - x_K|$.

The Lagrangian function associated with minimization problem (25) is naturally given by:

$$ \mathcal{L}_{\hat{K} \to \hat{N}}(U, \hat{U}, \lambda_1, \lambda_2) = \frac{1}{2} \sum_{K \in \hat{T}} m_K \left( u_K - \Pi_{\hat{K},K}(\hat{U})(x_K) \right)^2 - \lambda_1 \left( \sum_{K \in \hat{T}} m_K u_K - \sum_{\hat{K} \in \hat{T}} m_{\hat{K}} \hat{u}_{\hat{K}} \right) \\ + \lambda_2 \left( \sum_{K \in \hat{T}} m_K u_K^2 - \sum_{\hat{K} \in \hat{T}} m_{\hat{K}} \hat{u}_{\hat{K}}^2 \right), \quad (29) $$
while the Karush-Kuhn-Tucker first order optimality conditions correspond to:

$$ \left| \begin{array}{l} (1 + 2\lambda_2)u_K = \Pi_{\mathcal{N},K}(\hat{U})(\boldsymbol{x}_K) + \lambda_1, \\ \sum_{K \in \mathcal{T}} m_K u_K = \sum_{\hat{K} \in \hat{\mathcal{T}}} m_{\hat{K}} \hat{u}_{\hat{K}}, \\ \sum_{K \in \mathcal{T}} m_K u_K^2 = \sum_{\hat{K} \in \hat{\mathcal{T}}} m_{\hat{K}} \hat{u}_{\hat{K}}^2. \end{array} \right. \qquad (30) $$

Fortunately, system (30) can be directly solved. Indeed, denoting $\mu = 1 + 2\lambda_2$ and injecting the first relation into the discrete conservation constraint, we get:

$$ S_{\mathcal{N} \to \mathcal{N}}^{1}(\hat{U}) + |\Omega|\lambda_1 = \mu S_{\mathcal{N}}^{1}(\hat{U}), \qquad (31) $$

where we have denoted:

$$ S_{\mathcal{N} \to \mathcal{N}}^{1}(\hat{U}) = \sum_{K \in \mathcal{T}} m_K \Pi_{\mathcal{N},K}(\hat{U})(\boldsymbol{x}_K) \quad \text{and} \quad S_{\mathcal{N}}^{1}(\hat{U}) = \sum_{\hat{K} \in \hat{\mathcal{T}}} m_{\hat{K}} \hat{u}_{\hat{K}}. $$

Now, injecting the expression of $u_K$ in the $L^2$ norm constraint and naturally denoting:

$$ S_{\mathcal{N} \to \mathcal{N}}^{2}(\hat{U}) = \sum_{K \in \mathcal{T}} m_K \left| \Pi_{\mathcal{N},K}(\hat{U})(\boldsymbol{x}_K) \right|^2 \quad \text{and} \quad S_{\mathcal{N}}^{2}(\hat{U}) = \sum_{\hat{K} \in \hat{\mathcal{T}}} m_{\hat{K}} \hat{u}_{\hat{K}}^2, $$

we get:

$$ S_{\mathcal{N} \to \mathcal{N}}^{2}(\hat{U}) + 2\lambda_1 S_{\mathcal{N} \to \mathcal{N}}^{1}(\hat{U}) + |\Omega|\lambda_1^2 = \mu^2 S_{\mathcal{N}}^{2}(\hat{U}). $$

Combined with (31), this rewrites:

$$ S_{\mathcal{N} \to \mathcal{N}}^{2}(\hat{U}) + \frac{2}{|\Omega|} \left( \mu S_{\mathcal{N}}^{1}(\hat{U}) - S_{\mathcal{N} \to \mathcal{N}}^{1}(\hat{U}) \right) S_{\mathcal{N} \to \mathcal{N}}^{1}(\hat{U}) \\ + \frac{1}{|\Omega|} \left( \mu^2 S_{\mathcal{N}}^{1}(\hat{U})^2 - 2\mu S_{\mathcal{N}}^{1}(\hat{U}) S_{\mathcal{N} \to \mathcal{N}}^{1}(\hat{U}) + S_{\mathcal{N} \to \mathcal{N}}^{1}(\hat{U})^2 \right) = \mu^2 S_{\mathcal{N}}^{2}(\hat{U}), $$

ultimately leading to the second order equation:

$$ \mu^2 \left( S_{\mathcal{N}}^{2}(\hat{U}) - \frac{1}{|\Omega|} S_{\mathcal{N}}^{1}(\hat{U})^2 \right) = S_{\mathcal{N} \to \mathcal{N}}^{2}(\hat{U}) - \frac{1}{|\Omega|} S_{\mathcal{N} \to \mathcal{N}}^{1}(\hat{U})^2. \qquad (32) $$

Using Cauchy-Schwarz inequality, notice that:

$$ S_{\mathcal{N} \to \mathcal{N}}^{1}(\hat{U})^2 \le |\Omega| S_{\mathcal{N} \to \mathcal{N}}^{2}(\hat{U}) \quad \text{and} \quad S_{\mathcal{N}}^{1}(\hat{U})^2 \le |\Omega| S_{\mathcal{N}}^{2}(\hat{U}). $$

In the special case where $|\Omega|S_{\mathcal{N}}^{2}(\hat{U}) = S_{\mathcal{N}}^{1}(\hat{U})^2$, we deduce from the fact that Cauchy-Schwarz inequality is in this case an equality that $\hat{U}$ is a constant vector and we can simply use this constant value to define $P_{\mathcal{N} \to \mathcal{N}}(\hat{U})$. Thus, apart from the very specific case where $|\Omega|S_{\mathcal{N}}^{2}(\hat{U}) = S_{\mathcal{N}}^{1}(\hat{U})^2$, we see that (32) has always two real opposite solutions and that $P_{\mathcal{N} \to \mathcal{N}}(\hat{U})$ can correctly be defined, for instance always choosing the positive solution. To summarize this, we can write:

$$ P_{\mathcal{N} \to \mathcal{N}}(\hat{U})_K = \begin{cases} \frac{1}{\text{card}(\hat{\mathcal{T}})} \sum_{\hat{K} \in \hat{\mathcal{T}}} \hat{u}_{\hat{K}} & \text{if } |\Omega|S_{\mathcal{N}}^{2}(\hat{U}) = S_{\mathcal{N}}^{1}(\hat{U})^2, \\ \frac{1}{\mu(\hat{U})} \left( \Pi_{\mathcal{N},K}(\hat{U})(\boldsymbol{x}_K) + \frac{1}{|\Omega|} \mu(\hat{U}) S_{\mathcal{N}}^{1}(\hat{U}) - \frac{1}{|\Omega|} S_{\mathcal{N} \to \mathcal{N}}^{1}(\hat{U}) \right) & \text{Otherwise,} \end{cases} \qquad (33) $$
with of course:

$$ \mu(\hat{U}) = \left( \frac{|\Omega| S_{\hat{N} \to \mathcal{N}}^2(\hat{U}) - S_{\hat{N} \to \mathcal{N}}^1(\hat{U})^2}{|\Omega| S_{\hat{N}}^2(\hat{U}) - S_{\hat{N}}^1(\hat{U})^2} \right)^{1/2} . \quad (34) $$

Clearly, the projection operator $P_{\hat{N} \to \mathcal{N}}$ is extremely easy to compute and can be applied with almost zero cost regarding the other steps of the method.

# 6. NUMERICAL EXPLORATION

In our numerical experiments, the stabilization matrices are given by:

$$ S_{K,\sigma,\sigma'} = \frac{1}{|\eta_{K,\sigma}|^2} \Lambda_K \eta_{K,\sigma} \cdot \eta_{K,\sigma} \delta_{\sigma,\sigma'} \quad \text{and} \quad S_{K,\sigma,\sigma'}^0 = \gamma \beta_K \delta_{\sigma,\sigma'} $$

where $\delta_{\sigma,\sigma'}$ is the Kronecker symbol and $\gamma$ a parameter taken as $\gamma = 0.01$ in the following experiments. The discretization networks are generated following the procedure described in [14] that is based on the node generator of [31], while the geometries are generated by solving the linear system for exactly or non exactly compatible geometries presented in [14] (see also Appendix A).

## 6.1. Validation of the conservative NEM on reference test cases

Inspired by [16], on the square domain $\Omega =]0, 1]^2$ and the cubic domain $\Omega =]0, 1]^3$, we solve for the analytic solutions:

$$ u(x, y) = \sin(\pi x)\sin(\pi y) \text{ in } 2d \quad \text{and} \quad u(x, y, z) = \sin(\pi x)\sin(\pi y)\sin(\pi z) \text{ in } 3d, $$

our advection-diffusion-reaction problem for the parameters $\Lambda = \kappa \mathbb{I}_d, \beta = 1$ and:

$$ \alpha = \left( \frac{1}{4} + \frac{1}{8}x^2 - \frac{1}{4}y, \frac{1}{4} + \frac{1}{8}y^2 - \frac{1}{4}x \right)^T \quad \text{in } 2d, $$

and

$$ \alpha = \left( \frac{1}{4} + \frac{1}{8}x^2 + \frac{1}{4}(y-z), \frac{1}{4} + \frac{1}{8}y^2 + \frac{1}{4}(z-x), \frac{1}{4} + \frac{1}{8}z^2 + \frac{1}{4}(x-z) \right)^T \quad \text{in } 3d, $$

with $\kappa$ taking values $\{0.001, 1, 1000\}$ corresponding to the advection-dominated, intermediate and diffusion dominated regimes. The main difference with the test cases of [16] is the quadratic terms in the velocity fields $\alpha$, which we consider here to get a non constant divergence $\operatorname{div} \alpha$.

In dimension 2, the cell points generated using the node generator of [31] are presented on Figure 2. To assess the robustness of the method with respect to the distortion of the point cloud, in dimension 2 we apply the following distortion introduced in [34] and already used in [31] for meshes to both our cell and interfaces nodes (except boundary interfaces):

$$ \begin{aligned} x' &= x + 7(1 - x - y)xy(1 - x)(1 - y), \\ y' &= y + 7(1 - x - y)xy(1 - x)(1 - y), \end{aligned} \quad (35) $$

leading to the distorted point cloud also presented on Figure 2. We see that the effect of transformation (35) is to concentrate nodes on the surface $x+y=1$. In dimension 3, the cell points generated using the node generator of [31] are presented in Figure 3, on which we see that the node distribution is already quite random.

To study the $L^2$ convergence of solutions, we will consider

$$ \left( \sum_{K \in \mathcal{T}} m_K |u(\bar{x}_K) - \mathcal{M}_K(U)|^2 \right)^{1/2}, $$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Cell nodes in dimension 2: left reference cloud, right distorted cloud (reproduced from [31]).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Cell nodes in dimension 3.</figcaption></figure>

which will be in principle a $h^2$ approximation of the true $L^2$ error, provided the exact solution $u$ is regular enough. In the same way, we will use

$$ \left( \sum_{K \in \mathcal{T}} m_K |\nabla u(\bar{x}_K) - \nabla_K(U)|^2 \right)^{1/2} \quad \text{and} \quad \left( \sum_{K \in \mathcal{T}} m_K |\alpha \cdot \nabla u(\bar{x}_K) - G_K(\alpha_{\mathcal{F}_K}, U_K)|^2 \right)^{1/2}, $$

to the study the $L^2$ convergence of the gradient and the advective derivative. For this validation test case, we use exactly compatible geometries to construct our conservative NEM scheme. On Figures 4–6 we display convergence curves with respect to $h$ for respectively the $L^2$ norm of the error, the error on the gradient and on the error on the advective derivative $\alpha \cdot \nabla$. The approximate orders of convergence computed from those curves are gathered in Table 1.

We see that for both $\kappa = 1$ and $\kappa = 1000$, the convergence rates in $L^2$ norm of the solution are close to 2, as was observed for pure-diffusion reaction in [13]. For $\kappa = 0.001$, we observe the same reduction of convergence
<figure><img src="image_5.png" /><figcaption>FIGURE 4. $L^2$ convergence curves for advection-diffusion-reaction in dimension 2 and 3.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. $L^2$ convergence curves of the gradient for advection-diffusion-reaction in dimension 2 and 3.</figcaption></figure>

order that was described in [16] for the mesh-based hybrid-high order scheme. For the 2d distorted networks, we see that the main effect of distortion is to increase the error constant. Notice that the apparent convergence order increases in Table 1 precisely because of the distortion: most cells of the network have a diameter $h_K$ much smaller than $h$, thus the heterogeneity of $h_K$ induces a higher apparent order. The situation is similar for the errors on the gradient and on the advective derivative: for $\kappa = 1$ and $\kappa = 1000$, we obtain orders of convergence slightly greater than one, while in the advection dominated case the convergence rates decrease to around one half. This is again coherent with the results obtained for mesh based methods. Nevertheless notice that if the convergence rate are quite similar for the gradient and the advective derivative, the magnitude of the
<figure><img src="image_7.png" /><figcaption>FIGURE 6. $L^2$ convergence curves of the advective derivative for advection-diffusion-reaction in dimension 2 and 3.</figcaption></figure>

<table><caption>TABLE 1. Approximate orders of convergence for advection-diffusion-reaction.</caption><thead><tr><th rowspan="2">κ</th><th colspan="3">||·||<sub>L<sup>2</sup></sub></th><th colspan="3">||∇·||<sub>L<sup>2</sup></sub></th><th colspan="3">||(α·∇)·||<sub>L<sup>2</sup></sub></th></tr><tr><th>0.001</th><th>1</th><th>1000</th><th>0.001</th><th>1</th><th>1000</th><th>0.001</th><th>1</th><th>1000</th></tr></thead><tbody><tr><td>2d</td><td>1.31</td><td>2.14</td><td>2.3</td><td>0.56</td><td>1.18</td><td>1.18</td><td>0.57</td><td>1.25</td><td>1.25</td></tr><tr><td>2d-distorted</td><td>1.22</td><td>2.8</td><td>2.7</td><td>0.49</td><td>1.29</td><td>1.32</td><td>0.49</td><td>1.44</td><td>1.48</td></tr><tr><td>3d</td><td>1.12</td><td>2.24</td><td>2.34</td><td>0.59</td><td>1.33</td><td>1.33</td><td>0.67</td><td>1.38</td><td>1.37</td></tr></tbody></table>

error is one order of magnitude greater than the gradient than for the advective derivative. This also strengthens the idea that vanishing diffusion robust error estimates could be obtained for the NEM and that they would provide the same order of convergence as their equivalent for mesh-based methods.

## 6.2. Lagrangian vs. Eulerian advection

To illustrate the potential of the conservative NEM to easily handle both Lagrangian formulations and moving domains, let us consider a time dependent advection-diffusion problem:

$$ \frac{\partial u}{\partial t}(\mathbf{x}, t) - \kappa \Delta u(\mathbf{x}, t) + \operatorname{div}(\alpha(\mathbf{x}, t)u(\mathbf{x}, t)) = f(\mathbf{x}, t) \quad \text{in } \Omega(t) \times ]0, T[, \qquad (36) $$

with inhomogeneous Dirichlet boundary conditions $u(\mathbf{x}, t) = g(\mathbf{x}, t)$ on $\partial\Omega(t)$, a constant diagonal tensor $\Lambda = \kappa \mathbb{I}_d$ and a domain $\Omega(t)$ that is the result of the motion of an initial domain $\Omega_0$ by an incompressible speed field $\alpha$ (i.e. $\operatorname{div}(\alpha) = 0$). We denote $\mathbf{X}(t, \cdot)$ the mapping describing the motion of $\Omega$, i.e. for any $\mathbf{x}_0 \in \Omega_0$, the position at time $t$ of the point of domain $\Omega$ originally at $\mathbf{x}_0$ is given by $\mathbf{X}(t, \mathbf{x}_0)$. Denoting $U(\mathbf{x}_0, t) = u(\mathbf{X}(t, \mathbf{x}_0), t)$ and $F(\mathbf{x}_0, t) = f(\mathbf{X}(t, \mathbf{x}_0), t)$ the associated Lagrangian formulation is of course given by:

$$ \frac{\partial U}{\partial t}(\mathbf{x}_0, t) - \kappa \Delta U(\mathbf{x}_0, t) = F(\mathbf{x}_0, t) \quad \text{in } \Omega_0 \times ]0, T[, \qquad (37) $$
which rewrites using the usual Lagrangian formalism

$$ \frac{du}{dt}(x,t) - \kappa \Delta u(x,t) = f(x,t) \quad \text{in } \Omega(t) \times ]0,T[ \quad \text{where} \quad \frac{du}{dt} = \frac{\partial u}{\partial t} + \boldsymbol{\alpha} \cdot \nabla u. \qquad (38) $$

Using an implicit time scheme and given a discretization of $]0, T[$ in $N_T$ subintervals $]t_n, t_{n+1}$ [with $t_0 = 0$, $t_{N_T} = T$ and $\Delta t_n = t_{n+1} - t_n$, the time semi discrete scheme associated with the eulerian formulation is, denoting $\Omega_n = \Omega(t^n)$:

$$ \frac{u^{n+1} - u^n}{\Delta t_n} - \kappa \Delta u^{n+1} + \text{div}(\boldsymbol{\alpha}^{n+1} u^{n+1}) = f^{n+1} \quad \text{in } \Omega_{n+1}, \qquad (39) $$

while the time semi-discrete scheme associated with the Lagrangian formulation is

$$ \left| \begin{array}{l} \frac{u^{n+1}(X^{n+1}) - u^n(X^n)}{\Delta t_n} - \kappa \Delta u^{n+1}(X^{n+1}) = f^{n+1}(X^{n+1}) \quad \text{in } \Omega_{n+1}, \\ \\ \frac{X^{n+1} - X^n}{\Delta t_n} = \boldsymbol{\alpha}(X^n, t^n) \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text{in } \Omega_{n+1}. \end{array} \right. \qquad (40) $$

Clearly, those semi-discrete problems are advection-diffusion-reaction problems on which we can apply our conservative NEM scheme. Taking advantage of the meshless nature of the conservative NEM, the Lagrangian formulation will simply require to define a network on the initial geometry $\Omega_0$ and then move the points of the network using the Lagrangian speed field $\boldsymbol{\alpha}$ as specified in the second equation of (40). We will nevertheless need to recompute a network geometry at each new time step, after having displaced the domain. In the Eulerian case we need to compute a new network and its associated geometry on each new geometry $\Omega_n$ at each time step, and project the solution defined on the old network onto the newly constructed one.

In the Eulerian setting, let us denote $(\mathcal{N}^n, \mathcal{G}^n)$ the discretization network and geometry associated with $\Omega_n$. Both need to be recomputed for each time step in this context. Discarding the reaction stabilization (and thus the reaction flux, which is equivalent to set $\gamma = 0$) for this time dependent case, our Eulerian conservative NEM scheme on the moving domain $\Omega(t)$ is given for the time interval $]t_n, t_{n+1}$ [associated with the network $\mathcal{N}^{n+1}$ and geometry $\mathcal{G}^{n+1}$ by:

$$ \left| \begin{array}{l} \frac{m_K^{n+1}}{\Delta t_n} (u_K^{n+1} - P_{\mathcal{N}^n \to \mathcal{N}^{n+1}}(\mathbf{U}^n)_K) + m_K^{n+1} \delta_K(\boldsymbol{\alpha}) u_K^{n+1} + \\ \qquad \sum_{\sigma \in \mathcal{T}_K^n} F_{K,\sigma}^{D,n+1}(\mathbf{U}_K^{n+1}) + F_{K,\sigma}^{A,n+1}(\mathbf{U}_K^{n+1}) = m_K^{n+1} f_K^{n+1} \quad \forall K \in \mathcal{T}^{n+1}, \\ \\ \sum_{K \in \mathcal{T}_\sigma^n} F_{\sigma,K}^{D,n+1}(\mathbf{U}_K^{n+1}) + F_{\sigma,K}^{A,n+1}(\mathbf{U}_K^{n+1}) = -\frac{1}{2} \boldsymbol{\alpha}_\sigma^{n+1} \cdot \left( \sum_{K \in \mathcal{T}_\sigma^n} \boldsymbol{\eta}_{K,\sigma}^{n+1} \right) u_\sigma^{n+1} \quad \forall \sigma \in \mathcal{F}_{\text{int}}^{n+1}, \\ \\ u_\sigma^{n+1} = g(\mathbf{x}_\sigma, t^{n+1}) \quad \forall \sigma \in \mathcal{F}_{\text{ext}}^{n+1}. \end{array} \right. \qquad (41) $$

In the Lagrangian setting let us simply denote $\mathcal{N}$ the network associated with $\Omega_0$. At each time step, the position of the points of the network are updated through the Lagrangian formula:

$$ \left| \begin{array}{l} \mathbf{x}_K^{n+1} = \mathbf{x}_K^n + \Delta t_n \boldsymbol{\alpha}(\mathbf{x}_K^n, t^n) \quad \forall K \in \mathcal{T}, \\ \\ \mathbf{x}_\sigma^{n+1} = \mathbf{x}_\sigma^n + \Delta t_n \boldsymbol{\alpha}(\mathbf{x}_\sigma^n, t^n) \quad \forall \sigma \in \mathcal{F}, \end{array} \right. \qquad (42) $$

but the connectivity of the network is left unchanged. This implies that we may lose the admissibility of the network if the speed field induces too many distortion. A geometry $\mathcal{G}^{n+1}$ needs to be recomputed for each time
new time step, leading to new measures $(m_K^{n+1})_{K \in \mathcal{T}}$ and new $(\eta_{K,\sigma}^{n+1})_{\sigma \in \mathcal{F}_K, K \in \mathcal{T}}$. Again discarding the reaction flux, the Lagrangian conservative NEM scheme on the moving domain $\Omega(t)$ is given for the time interval $]t_n, t_{n+1}$[ by, using the moving network $\mathcal{N}$ whose geometry for this time interval is $\mathcal{G}^{n+1}$:

$$ 
\left\{
\begin{array}{l}
\displaystyle \frac{m_K^{n+1}}{\Delta t_n} (u_K^{n+1} - u_K^n) + \sum_{\sigma \in \mathcal{F}_K} F_{K,\sigma}^{D,n+1}(\mathbf{U}_K^{n+1}) = m_K^{n+1} f_K^{n+1} \quad \forall K \in \mathcal{T}, \\[1.5em]
\displaystyle \sum_{K \in \mathcal{T}_\sigma} F_{\sigma,K}^{D,n+1}(\mathbf{U}_K^{n+1}) = 0 \qquad \forall \sigma \in \mathcal{F}_{\text{int}}, \\[1.5em]
\displaystyle u_\sigma^{n+1} = g(\mathbf{x}_\sigma, t^{n+1}) \qquad \forall \sigma \in \mathcal{F}_{\text{ext}}.
\end{array}
\right.
\qquad (43)
$$

**Algorithm 2 Eulerian advection**

1: Initialization:

2: Compute $(\mathcal{N}^0, \mathcal{G}^0)$ for domain $\Omega_0$.

3: Compute $\mathbf{U}^0$ on $\mathcal{N}^0$.

4: Time loop:

5: **for** $n = 0...N_T - 1$ **do**

6: Move $\Omega_n$ into $\Omega_{n+1}$.

7: Compute $(\mathcal{N}^{n+1}, \mathcal{G}^{n+1})$ for domain $\Omega_{n+1}$.

8: Compute the projection $P_{\mathcal{N}^{n-1}\rightarrow\mathcal{N}^{n+1}}(\mathbf{U}^n)$ of $\mathbf{U}^n$ onto $\mathcal{N}^{n+1}$.

9: Solve (41) to compute $\mathbf{U}^{n+1}$.

10: **end for**

**Algorithm 3 Lagrangian advection**

1: Initialization:

2: Compute $(\mathcal{N}, \mathcal{G}^0)$ for domain $\Omega_0$.

3: Compute $\mathbf{U}^0$ on $\mathcal{N}$.

4: Time loop:

5: **for** $n = 0...N_T - 1$ **do**

6: Move $\Omega_n$ into $\Omega_{n+1}$.

7: Move $\mathcal{N}$ using (42).

8: Compute $\mathcal{G}^{n+1}$ for domain $\Omega_{n+1}$.

9: Solve (43) to compute $\mathbf{U}^{n+1}$.

10: **end for**

In Algorithms 2 and 3, we respectively summarize the needed steps for the Eulerian and Lagrangian approaches. We see that the network geometry has to be recomputed at each time step in both cases, the Lagrangian approach only saving the network generation. However the advection is handled in the Lagrangian case directly through the motion of the network. In the NEM case the main advantage of the Lagrangian approach is the natural handling of advection and, if one prefers to use explicit schemes for the advection, no potentially stringent advective CFL constraint.

Now, we consider a slight modification of the 2d test case of [30]. The domain $\Omega_0 = B_1 \setminus (E_1 \cup E_2)$ corresponds to the unit ball minus the two sets $E_1$ and $E_2$ given in parametric form by:

$$ 
\begin{aligned}
E_1: & \quad x = 0.4 \cos(\theta) \quad y = -0.5 + 0.2 \sin(\theta) \quad \theta \in [0, 2\pi[, \\
E_2: & \quad x = 0.1 \cos(\theta) \quad y = 0.3 + 0.2 \sin(\theta) \quad \theta \in [0, 2\pi[.
\end{aligned}
$$

The incompressible speed field is given by

$$ \alpha_x(\mathbf{x}, t) = 2 \sin(\pi |\mathbf{x}|^2) \sin(\pi t)y \quad \alpha_y(\mathbf{x}, t) = -2 \sin(\pi |\mathbf{x}|^2) \sin(\pi t)x, \qquad (44) $$

while we seek the exact solution given by:

$$ u(\mathbf{x}, t) = (1 + \sin(\pi x) \cos(\pi y) \sin(\pi t)). \qquad (45) $$

Notice that the speed field cancels on the boundary of $B_1$, which will increase the distortion near $\partial B_1$. We consider the two values $\kappa = 1$ and $\kappa = 0.001$ to again explore the behavior up to the advection dominated case. We solve from $t_0 = 0$ to $T = 1$ with a constant time step $\Delta t=0.01$, transforming the originally elliptic holes into the quite deformed ones displayed on Figure 7.

Contrary to the Eulerian case, to enforce exact conservation in the Lagrangian case as we formally get a diffusion-reaction system we have seen that there is no need for exactly compatible geometries. Consequently
<figure><img src="image_8.png" /><figcaption>FIGURE 7. The moving domain with holes. Left: initial state at $t = 0$. Right: final state at $t = 1$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 8. Relative $L^2$ error with respect to time for time dependent advection-diffusion on the moving domain with holes.</figcaption></figure>

we do not enforce exact compatibility in the Lagrangian geometry computation. On Figure 8, we display the evolution with time of the relative $L^2$ error (i.e. the $L^2$ error divided by the $L^2$ norm of the exact solution). Clearly, in the advection dominated case ($\kappa = 0.001$) the Lagrangian approach produces better results, while in the diffusion dominated case ($\kappa = 1$) the Lagrangian approach is much less satisfactory. This is exactly the opposite for the Eulerian approach: in the advection dominated case it leads to larger errors than the Lagrangian approach and the contrary in the diffusion dominated case. Our understanding of those results is quite classical: in the advection-dominated case, the numerical diffusion induced by the projection $P_{N^n \rightarrow N^{n+1}}$ accumulates over time for the Eulerian approach, finally generating too many diffusion, while in the diffusion dominated case the numerical diffusion remains negligible when compared to the true diffusion. Thus, quite naturally the Lagrangian approach which does not suffer from numerical diffusion produces much better results
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Lagrangian evolution of the network for the time dependent advection-diffusion case. *Left*: initial network, *right*: final network.</figcaption></figure>

in the advection dominated case. In the diffusion dominated case, the origin of the large error of the Lagrangian approach in the diffusion dominated case becomes obvious when looking at the evolution of the network with time presented on Figure 9: clearly, the network obtained by a simple Lagrangian displacement of the cell and interface nodes becomes too distorted to get a geometry with enough quality to correctly represent the gradients and the diffusion operator. The error peak observed in Figure 8 coincides with the moment where the network reaches this high level of distortion, the decrease observed afterwards being due to the fact that the exact solution goes back to a constant and is thus easier to approximate even on a distorted network. Looking at the final network of Figure 9, we see that most of the distortion comes from the interfaces movement. As moreover scheme (43) only transports cell values from one time step to another and as we have to recompute the geometry anyway, a simple remedy would be to redefine only the interfaces at each time step before computing the geometry. This corresponds to the “Lagrange+” results of Figure 8. We see that this clearly solves the distortion issue in this particular case, confirming that the large error of the pure Lagrangian approach did come from network distortion.

We conclude our study of this time dependent case by displaying some cumulative computational times in Table 2. Obviously, as we only generate the network once in the pure Lagrangian case, the cost of network generation is completely negligible in this case. Comparing the cost of network generation in the Eulerian and “Lagrange+” case, we observe that interface generation is the most costly part of the algorithm proposed in [14]. We retrieve the fact explained in [14] that the main part to improve in the network generation procedure of is the interface generation, and not the cell generation that is based on the node generator of [31]. We also see on Table 2 that the most costly part of the overall NEM workflow is the geometry generation. As it simply corresponds to solving a linear system, this implies that our combination of CG with SSOR preconditioning while very robust could certainly be replaced by a more efficient dedicated solver. This also strongly advocates for the use of the NEM for larger non-linear systems for which the solving time will be prominent. Notice that we have used our own implementation of CG and SSOR in the prototype code that supports all those numerical experiments. As it is far from being as highly optimized as can be alternatives coming from state of the art linear algebra libraries, our results in this regard are necessarily quite pessimistic.
<table><caption>TABLE 2. Cumulative computational times (in seconds) for the time dependent case with moving holes.</caption><thead><tr><th rowspan="2">κ</th><th colspan="2">Network generation</th><th colspan="2">Geometry generation</th><th colspan="2">Solver</th></tr><tr><th>0.001</th><th>1</th><th>0.001</th><th>1</th><th>0.001</th><th>1</th></tr></thead><tbody><tr><td>Euler</td><td>9.327</td><td>9.248</td><td>87.829</td><td>86.833</td><td>3.457</td><td>3.797</td></tr><tr><td>Lagrange</td><td>0.21</td><td>0.128</td><td>64.78</td><td>63.78</td><td>3.9150</td><td>3.68</td></tr><tr><td>Lagrange+</td><td>7.32</td><td>6.94</td><td>95.48</td><td>97.51</td><td>4.55</td><td>5.16</td></tr></tbody></table>

### 6.3. Application to a simple two-phase Darcy flow model

In this subsection we try to illustrate the interest of the NEM for basin simulations, through the NEM Lagrangian simulation of a simplified two-phase Darcy flow model on a moving basin. Of course, as our aim here is purely to give an idea of the potential of the NEM in such situations, we will consider a relatively simple yet realistic basin geometry, that could obviously be meshed with relative ease. Generating networks for complex faulted basin indeed goes far beyond the scope of the present study. Nevertheless, this will also allow us to explore for the first time the behavior of the conservative NEM (and more generally of the NEM) in a non-linear setting.

Let us consider the following simplified two-phase Darcy flow model:

$$ 
\left\{
\begin{array}{l}
\phi \frac{\partial(1-s)}{\partial t} - \operatorname{div} \left( \frac{\rho_w}{\mu_w} k_{r,w}(s) \Lambda(\nabla P_w - \rho_w g \zeta) \right) = f_w \quad \text{in } \Omega, \\
\\
\phi \frac{\partial s}{\partial t} - \operatorname{div} \left( \frac{\rho_{nw}}{\mu_{nw}} k_{r,nw}(s) \Lambda(\nabla P_{nw} - \rho_{nw} g \zeta) \right) = f_{nw} \quad \text{in } \Omega, \\
\\
P_{nw} = P_w + p_c(s) \quad \text{in } \Omega,
\end{array}
\right.
\quad (46) $$

with boundary conditions:

$$ 
\left\{
\begin{array}{ll}
s = s_b & \text{on } \partial\Omega_D \\
\\
P_w = P_b & \text{on } \partial\Omega_D \\
\\
\frac{\rho_w}{\mu_w} k_{r,w}(s) \Lambda(\nabla P_w - \rho_w g \zeta) \cdot \mathbf{n} = 0 & \text{on } \partial\Omega_N, \\
\\
\frac{\rho_{nw}}{\mu_{nw}} k_{r,nw}(s) \Lambda(\nabla P_{nw} - \rho_{nw} g \zeta) \cdot \mathbf{n} = 0 & \text{on } \partial\Omega_N,
\end{array}
\right.
\quad (47) $$

where we have split the boundary of $\Omega$ in two parts $\partial\Omega = \overline{\partial\Omega_D} \cup \overline{\partial\Omega_N}$, with $\partial\Omega_D$ corresponding to Dirichlet boundary conditions and $\partial\Omega_N$ corresponding to Neumann ones, where $\zeta$ is the depth, $\Lambda$ is the permeability tensor of the porous media, $\phi$ its porosity, $s$ represents the saturation of a non wetting phase, $P_w$ and $P_{nw}$ the pressures of respectively the wetting and non-wetting phases, $k_{r,w}$ the relative permeability law of the wetting phase, $k_{r,nw}$ and $p_c$ the relative permeability and capillary pressure laws of the non-wetting phase, and $(\rho_w, \mu_w)$ and $(\rho_{nw}, \mu_{nw})$ the constant density and viscosity of the wetting and non wetting phases.
We assume here for simplicity that the discretization network $\mathcal{N}$ has only “face like” interfaces, i.e. for every $\sigma \in \mathcal{F}_{\text{int}}$, $\text{card}(\mathcal{T}_{\sigma}) = 2$ and for every $\sigma \in \mathcal{F}_{\text{ext}}$, $\text{card}(\mathcal{T}_{\sigma}) = 1$. This precisely corresponds to the output network of the procedure described in [14] that is based on the node generator of [31]. Then, we consider the following Lagrangian discretization using the conservative NEM scheme and a fully implicit time scheme:

$$ 
\begin{aligned}
& \frac{m_K^{n+1} \phi_K}{\Delta t_n} (s_K^{n+1} - s_K^n) + \sum_{\sigma \in \mathcal{F}_K \setminus \mathcal{F}_{\text{ext}}^N} \frac{\rho_{nw}}{\mu_{nw}} k_{r,nw,\sigma}^{n+1} F_{K,\sigma}^{D,n+1} (\mathbf{P}_{nw,K}^{n+1} - \rho_{nw} g \zeta_K) = m_K^{n+1} f_w(\mathbf{x}_K, t^{n+1}) \quad \forall K \in \mathcal{T}, \\
& \sum_{\sigma \in \mathcal{F}_K \setminus \mathcal{F}_{\text{ext}}^N} \frac{\rho_w}{\mu_w} k_{r,w,\sigma}^{n+1} F_{K,\sigma}^{D,n+1} (\mathbf{P}_{w,K}^{n+1} - \rho_w g \zeta_K) + \frac{\rho_{nw}}{\mu_{nw}} k_{r,nw,\sigma}^{n+1} F_{K,\sigma}^{D,n+1} (\mathbf{P}_{nw,K}^{n+1} - \rho_{nw} g \zeta_K) \\
& = m_K^{n+1} (f_w(\mathbf{x}_K, t^{n+1}) + f_{nw}(\mathbf{x}_K, t^{n+1})) \quad \forall K \in \mathcal{T}, \\
& P_{nw,K}^{n+1} = P_{w,K}^{n+1} + p_c(s_K^{n+1}) \quad \forall K \in \mathcal{T}, \\
& \sum_{K \in \mathcal{T}_{\sigma}} F_{\sigma,K}^{D,n+1} (\mathbf{P}_{w,K}^{n+1} - \rho_w g \zeta_K) = 0 \quad \forall \sigma \in \mathcal{F}_{\text{int}}, \\
& \sum_{K \in \mathcal{T}_{\sigma}} F_{\sigma,K}^{D,n+1} (\mathbf{P}_{nw,K}^{n+1} - \rho_{nw} g \zeta_K) = 0 \quad \forall \sigma \in \mathcal{F}_{\text{int}},
\end{aligned}
\quad (48) $$

along with the boundary conditions:

$$ 
\left\{
\begin{array}{ll}
s_{\sigma}^{n+1} = s_b(\mathbf{x}_{\sigma}, t^{n+1}) & \forall \sigma \in \mathcal{F}_{\text{ext}}^{D}, \\
P_{w,\sigma}^{n+1} = P_b(\mathbf{x}_{\sigma}, t^{n+1}) & \forall \sigma \in \mathcal{F}_{\text{ext}}^{D}, \\
P_{nw,\sigma}^{n+1} = P_{w,\sigma}^{n+1} + p_c(s_{\sigma}^{n+1}) & \forall \sigma \in \mathcal{F}_{\text{ext}}^{D}, \\
\sum_{K \in \mathcal{T}_{\sigma}} F_{\sigma,K}^{D,n+1} (\mathbf{P}_{w,K}^{n+1} - \rho_w g \zeta_K) = 0 & \forall \sigma \in \mathcal{F}_{\text{ext}}^{N}, \\
\sum_{K \in \mathcal{T}_{\sigma}} F_{\sigma,K}^{D,n+1} (\mathbf{P}_{nw,K}^{n+1} - \rho_{nw} g \zeta_K) = 0 & \forall \sigma \in \mathcal{F}_{\text{ext}}^{N}.
\end{array}
\right.
\quad (49) $$

where we assume that $\mathcal{F}_{\text{ext}} = \mathcal{F}_{\text{ext}}^{D} \cup \mathcal{F}_{\text{ext}}^{N}$ and $\mathbf{x}_{\sigma} \in \partial \Omega_{D} \quad \forall \sigma \in \mathcal{F}_{\text{ext}}^{D}$ and $\mathbf{x}_{\sigma} \in \partial \Omega_{N} \quad \forall \sigma \in \mathcal{F}_{\text{ext}}^{N}$. Notice that we enforce a “linear” continuity only of the phase pressures on interior interfaces for simplicity. The assumption on the network allows finally to upwind the relative permeabilities on each interface $\sigma$ according to each phase flux, i.e. if $\mathcal{T}_{\sigma} = (K, L)$:

$$ k_{r,nw,\sigma}^{n+1} = \begin{cases} k_{r,nw}(s_K^{n+1}) & \text{if } F_{K,\sigma}^{D}(\mathbf{P}_{nw,K}^{n+1} - \rho_{nw}g\zeta_K) \ge 0, \\ k_{r,nw}(s_L^{n+1}) & \text{if } F_{K,\sigma}^{D}(\mathbf{P}_{nw,K}^{n+1} - \rho_{nw}g\zeta_K) < 0, \end{cases} \quad (50) $$

and

$$ k_{r,w,\sigma}^{n+1} = \begin{cases} k_{r,nw}(s_K^{n+1}) & \text{if } F_{K,\sigma}^{D}(\mathbf{P}_{w,K}^{n+1} - \rho_w g \zeta_K) \ge 0, \\ k_{r,w}(s_L^{n+1}) & \text{if } F_{K,\sigma}^{D}(\mathbf{P}_{w,K}^{n+1} - \rho_w g \zeta_K) < 0. \end{cases} \quad (51) $$

We consider the moving basin displayed on Figure 10, that initially ranges from $x = 0$ km to $x = 4.9723$ km and $y = -0.27783$ km to $y = -3.6556$ km. Starting from horizontally deposited layers, the right part is uplifted from 0 to 1 million years (My) until two anticlinal areas are formed where the non wetting phase can accumulate.
<figure><img src="image_11.png" /><figcaption>FIGURE 10. The moving porous basin. Top: initial state at $t = 0$ My. Mid left: state at $t = 0.5$ My. Mid right: state at $t = 1.0$ My. Bottom left: state at $t = 1.5$ My. Bottom right: state at $t = 2$ My.</figcaption></figure>

Next, from 1 to 2 My, we mirror the uplift to form the two same accumulation zones but on the left part of the basin. We produce the non wetting phase in the window $x = 3$ km to $x = 4.5$ km and $y = -3$ km to $y = -3.6$ km, from $t = 0.2$ My to $t = 0.8$ My, with a production linearly growing with time from 0 at $t = 0.2$ My until reaching its maximum at $t = 0.5$ My and then linearly going back to zero with time at $t = 0.8$ My. Our boundary conditions are open top boundary (inhomogenous Dirichlet) with $s_b = 0$ and wetting pressure $P_b$ set to the hydrostatic pressure with reference 1 bar at $z = 0$, and close (homogeneous Neumann) for the three other parts of the boundary. The density and viscosity of the wetting and non wetting phases are taken constant with
<figure><img src="image_12.png" /><figcaption>FIGURE 11. Porosity and $log_{10}$ of permeability (in milli Darcy).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12. Saturations for different times. Upper left $t = 0.3$ My. Upper right $t = 0.4$ My. Bottom left $t = 0.5$ My. Bottom right $t = 0.8$ My.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIGURE 13. Saturations for different times. Upper left $t = 1.0$ My. Upper right $t = 1.4$ My. Bottom left $t = 1.5$ My. Bottom right $t = 2.0$ My.</figcaption></figure>

values:

$$ \rho_w = 1000 \text{ kg.m}^{-3}, \quad \rho_{nw} = 700 \text{ kg.m}^{-3}, \quad \mu_w = 0.001 \text{ m}^2 \cdot \text{s}^{-1}, \quad \mu_{nw} = 0.005 \text{ m}^2 \cdot \text{s}^{-1}. $$

The gravity constant is $g = 9.81 \text{ m.s}^{-2}$, while we illustrate on Figure 11 the constant values of the porosity and permeability on the domain at time $t = 1$ My. Finally, the relative permeability laws and capillary pressure laws are chosen as simplified Brooks-Corey laws:

$$ k_{r,w}(s) = 1 - s^2, \quad k_{r,nw} = s^2, \quad p_c(s) = p_{c,e} + s^2, $$

where the entry capillary pressure is chosen as equal to zero everywhere except in the two layers with the lowest permeability values illustrated on Figure 11, where we choose $p_{c,e} = 10$ MPa, to form a barrier that will trap the non wetting phase.

We show the non wetting saturation for different times on Figures 12 and 13. As we can see, the non wetting phase has the expected behavior: it first moves vertically to accumulate in the anticlinal below the barrier, until it becomes too much concentrated and then spills from the bottom left into the middle anticlinal. Then, as the basin evolves the non wetting phase migrates to the left to accumulate under the newly formed anticlinals. This qualitative study confirms the potential of the NEM to be cope with complex Lagrangian simulations. To conclude, we show the computational times for this test case in Table 3. As we can see, contrary to the scalar
<table><caption>TABLE 3. Cumulative computational times (in seconds) for the two phase Darcy flow case.</caption><thead><tr><td>Geometry generation</td><td>Non linear solver</td></tr></thead><tbody><tr><td>194.68</td><td>956.48</td></tr></tbody></table>

and linear time dependent case, this time the geometry generation only takes 20% of the non linear solving time.

## 7. CONCLUSION AND PERSPECTIVES

We explained how to derive a conservative network element method for advection-diffusion-reaction problems which admits a finite volume like flux formulation. Numerical results in dimension 2 and 3 illustrate the good behavior and robustness of the method in both the diffusion and advection dominated regimes, even on distorted or moving point clouds. As the method can be formulated as a finite volume method one of its main interesting features is that its extension to more complex and general problems such as multiphase flow in porous media has been shown to follow simple, well established lines.

## APPENDIX A. NETWORK GEOMETRY GENERATION

In this section, we briefly recall the geometry generation procedure described in [14]. Consequently, all this section is directly reproduced from [14]. Admissible network geometries are obtained by minimizing the cost function:

$$ \mathcal{J}(\mathcal{G}) = \frac{\delta_{\mathcal{J}}^m}{2} \sum_{K \in \mathcal{T}} m_K^2 + \frac{\delta_{\mathcal{J}}^{\eta}}{2} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} |\eta_{K,\sigma}|^2 + \frac{1}{2} \sum_{K \in \mathcal{T}} \omega_{\mathcal{J},K}^0 |\boldsymbol{\epsilon}_{m,K}^0|^2 + \frac{1}{2} \sum_{K \in \mathcal{T}} \sum_{i=1}^d \sum_{j=1}^d \omega_{\mathcal{J},K}^1 |\boldsymbol{\epsilon}_{m,K}^{1,ij}|^2 + \frac{1}{2} \sum_{\sigma \in \mathcal{F}_{\text{int}}} \omega_{\mathcal{J},\sigma} |\boldsymbol{\epsilon}_{\sigma}|^2, \quad (A.1) $$

where $\delta_{\mathcal{J}}^m > 0$ and $\delta_{\mathcal{J}}^{\eta} > 0$ are parameters, where $\boldsymbol{\epsilon}_{m,K}^0 = m_K \boldsymbol{\epsilon}_K^0$ and $\boldsymbol{\epsilon}_{m,K}^{1,ij} = m_K \boldsymbol{\epsilon}_K^{1,ij}$, and where the weights $\omega_{\mathcal{J}}$ are chosen as $\omega_{\mathcal{J},K}^0 = h_K^{-(\nu+1)}$, $\omega_{\mathcal{J},K}^1 = h_K^{-(\nu+2)}$ and $\omega_{\mathcal{J},\sigma} = (\min_{K \in \mathcal{T}_{\sigma}} h_K^{\nu+1})^{-1}$ with $\nu > 0$ is parameter. This minimization is performed on the constrained set $\mathcal{A}_{\text{geom}}$ enforcing the consistency and conservation properties of the geometry:

$$ \mathcal{A}_{\text{geom}} = \left\{ \mathcal{G} \text{ such that } \sum_{K \in \mathcal{T}} m_K = |\Omega|, \text{ and } \sum_{\sigma \in \mathcal{F}_K} \eta_{K,\sigma}^i = m_K \boldsymbol{\epsilon}_K^{0,i} \quad \text{for all } K \in \mathcal{T} \text{ and all } 1 \le i \le d, \right. \\ \text{and } \sum_{\sigma \in \mathcal{F}_K} \eta_{K,\sigma}^i (x_{\sigma}^j - x_K^j) = m_K (\delta_{ij} + \boldsymbol{\epsilon}_K^{1,ij}) \quad \text{for all } K \in \mathcal{T} \text{ and all } 1 \le i, j \le d, \\ \left. \text{and } \sum_{K \in \mathcal{T}_{\sigma}} \eta_{K,\sigma}^i = \boldsymbol{\epsilon}_{\sigma}^i \quad \text{for all } \sigma \in \mathcal{F}_{\text{int}} \text{ and all } 1 \le i \le d \right\}. $$

The Karush-Kuhn-Tucker first order optimality conditions associated with this minimization problem corresponds to a linear system from which many unknowns can be eliminated through Schur's complement. This finally leads to the decomposition into a primary and a secondary system depicted in Table A.1, the secondary system being the part handled through Schur's complement.

To get an exactly compatible geometry, it suffices to take $\boldsymbol{\epsilon}_{\sigma}^i = 0$ and to remove the equation in front of $\boldsymbol{\epsilon}_{\sigma}^i$ in the system of Table A.1.
<table><caption>TABLE A.1. Equations for geometry generation (reproduced from [14]).</caption><thead><tr><th></th><th>Equation</th><th>Unknown</th></tr></thead><tbody><tr><td rowspan="2">Primary system</td><td>Σ<sub>K∈T<sub>σ</sub></sub> η<sup>i</sup><sub>Kσ</sub> = ε<sup>i</sup><sub>σ</sub></td><td>λ<sup>i</sup><sub>σ</sub></td></tr><tr><td>Σ<sub>K∈T</sub> m<sub>K</sub> = |Ω|</td><td>λ<sub>Ω</sub></td></tr><tr><td rowspan="7">Secondary system</td><td>ω<sub>J,σ</sub>ε<sup>i</sup><sub>σ</sub> - λ<sup>i</sup><sub>σ</sub> = 0</td><td>ε<sup>i</sup><sub>σ</sub></td></tr><tr><td>ω<sup>0</sup><sub>J,K</sub>ε<sup>0,i</sup><sub>m,K</sub> - λ<sup>0,i</sup><sub>K</sub> = 0</td><td>λ<sup>0,i</sup><sub>K</sub></td></tr><tr><td>ω<sup>1</sup><sub>J,K</sub>ε<sup>1,ij</sup><sub>m,K</sub> - λ<sup>1,ij</sup><sub>K</sub> = 0</td><td>λ<sup>1,ij</sup><sub>K</sub></td></tr><tr><td>Σ<sub>σ∈F<sub>K</sub></sub> η<sup>i</sup><sub>K,σ</sub> = ε<sup>0,i</sup><sub>m,K</sub></td><td>ε<sup>0,i</sup><sub>m,K</sub></td></tr><tr><td>Σ<sub>σ∈F<sub>K</sub></sub> η<sup>i</sup><sub>K,σ</sub>(x<sup>j</sup><sub>σ</sub> - x<sup>j</sup><sub>K</sub>) = m<sub>K</sub>δ<sub>ij</sub> + ε<sup>1,ij</sup><sub>m,K</sub></td><td>ε<sup>1,ij</sup><sub>m,K</sub></td></tr><tr><td>δ<sup>η</sup><sub>J</sub>η<sup>i</sup><sub>K,σ</sub> + Σ<sub>j=1</sub><sup>d</sup> λ<sup>1,ij</sup><sub>K</sub>(x<sup>j</sup><sub>σ</sub> - x<sup>j</sup><sub>K</sub>) + λ<sup>0,i</sup><sub>K</sub> + λ<sup>i</sup><sub>σ</sub> = 0</td><td>η<sup>i</sup><sub>K,σ</sub></td></tr><tr><td>δ<sup>m</sup><sub>J</sub>m<sub>K</sub> - Σ<sub>i=1</sub><sup>d</sup> λ<sup>1,ii</sup><sub>K</sub> + λ<sub>Ω</sub> = 0</td><td>m<sub>K</sub></td></tr></tbody></table>

## APPENDIX B. ELEMENTS OF NETWORK ELEMENT CONVERGENCE THEORY

In this section, we recall the main ingredients of the NEM convergence theory relying on the coercivity of the diffusion operator that was developed in [11–13]. As for the previous appendix, most of the results presented in this appendix are reproduced from [11–13] for the reader’s convenience. The key idea to bridge the gap between the degrees of freedom space $X_N$ on which the NEM is defined and the usual function spaces is the notion of quadrature families of [11], which is a family of functions $(\psi_K)_{K \in T}$ such that for any $K \in T$, $\psi_K \in L^\infty(\mathbb{R}^d)$ and:

$$ \int_{\Omega} \psi_K = m_K, \quad \sum_{K \in T} \psi_K = 1 \quad \text{for a.e } \mathbf{x} \in \Omega, \quad \text{supp } \psi_K \subset B(\mathbf{x}_K, \rho_K) \supset B_K, $$

We denote $B_K = B(\mathbf{x}_K, \rho_K)$. Quadrature families are the cornerstone of the convergence theory developed in [11], where an existence result for them was established. True functions are then reconstructed on $\mathbb{R}^d$ and $\Omega$ from the dofs using the $\psi_K$'s through:

$$ \Pi_T(\mathbf{U}) = \sum_{K \in T} \psi_K \mathcal{M}_K(\mathbf{U}) \quad \text{and} \quad \nabla_T(\mathbf{U}) = \sum_{K \in T} \psi_K \nabla_K(\mathbf{U}) \quad \text{and} \quad \Pi_N(\mathbf{U}) = \sum_{K \in T} \psi_K \Pi_K(\mathbf{U}_K). $$

and

$$ (\alpha_{\mathcal{F}} \cdot \nabla_{T})(\mathbf{U}) = \sum_{K \in T} \psi_K \mathbf{G}_K(\alpha_{\mathcal{F}_K}, \mathbf{U}). $$

If $(\mathcal{N}, \mathcal{G})$ is an admissible discretization network and an associated admissible network geometry, and if $\Omega$ satisfies the cone condition, we know from [11] that there exists $C > 0$ depending only on $\Omega$ and the quality parameters of geometry and of the quadrature family such that for any $\mathbf{U} \in X_N$:

$$ ||\Pi_T(\mathbf{U})||_{L^2(\Omega)} \le C ||\mathbf{U}||_X \quad \text{and} \quad ||\nabla_T(\mathbf{U})||_{L^2(\Omega)} \le C ||\mathbf{U}||_X \quad \text{and} \quad ||(\alpha_{\mathcal{F}} \cdot \nabla_T)(\mathbf{U})||_{L^2(\Omega)} \le C ||\mathbf{U}||_X, $$
where the discrete $H^1$ norm $||\cdot||_X$ on $X_N$ is defined by:

$$ ||U||_X^2 = (U, U)_X \quad \text{and} \quad (U, V)_X = \sum_{K \in \mathcal{T}} (U, V)_{0,K} + \sum_{K \in \mathcal{T}} (U, V)_{1,K}, $$

where

$$ (U, V)_{0,K} = m_K \mathcal{M}_K(U_K) \mathcal{M}_K(V_K) + \sum_{\sigma \in \mathcal{F}_K} m_K (u_\sigma - \mathcal{M}_K(U_K))(v_\sigma - \mathcal{M}_K(V_K)), $$

and

$$ (U, V)_{1,K} = \sum_{\sigma \in \mathcal{F}_K} m_K h_K^{-2} (u_\sigma - \mathcal{M}_K(U_K))(v_\sigma - \mathcal{M}_K(V_K)). $$

Conversely, to go from function spaces to degrees of freedom spaces, we recall the interpolation procedure described in [11]. To any $\sigma \in \mathcal{F}$, we associate a radius $r_\sigma > 0$ such that $B_\sigma \subset B_K$, where $B_\sigma = B(x_\sigma, r_\sigma)$. As $\Omega$ is assumed Lipschitz, using Stein's extension theorem [32] we also know that there exists an operator $E$ such that for any $k \ge 0$, there exists $C_{E,k} > 0$ such that for any $v \in H^k(\Omega)$, $Ev \in H^k(\mathbb{R}^d)$, $Ev = v$ in $\Omega$ and

$$ |Ev|_{H^k(\mathbb{R}^d)} \le C_{E,k} |v|_{H^k(\Omega)}, $$

and if $v \in H_0^1(\Omega)$, then $Ev = 0$ in $\mathbb{R}^d \setminus \Omega$. Then, we define the operator $\mathcal{I} : H^1(\Omega) \to X_N$ by setting $\mathcal{I}(v) = (\mathcal{I}_\sigma(v))_{\sigma \in \mathcal{F}}$ where:

$$ \mathcal{I}_{\sigma}(v) = \frac{1}{|B_{\sigma}|} \int_{B_{\sigma}} Ev \quad \text{for any } \sigma \in \mathcal{F}, $$

and we of course denote $\mathcal{I}_K(v) = (\mathcal{I}_\sigma(v))_{\sigma \in \mathcal{F}_K}$. We also introduce the operator $\mathcal{I}^0 : H_0^1(\Omega) \to X_N$

$$ \mathcal{I}_{\sigma}^{0}(v) = \begin{cases} \frac{1}{|B_{\sigma}|} \int_{B_{\sigma}} Ev & \text{for any } \sigma \in \mathcal{F}_{\text{int}} \\ 0 & \text{for any } \sigma \in \mathcal{F}_{\text{ext}}. \end{cases} $$

We will now briefly sketch the proof of an abstract error estimate. The first step is to notice that the stabilization terms are strongly consistent as one easily shows that (adapting the proof of lemma 5.11 of [12]), for any $v \in H^2(\Omega)$:

$$ s_h(\mathcal{I}(v), \mathcal{I}(v)) \le CS^* \left( \sum_{K \in \mathcal{T}} h_K^{-2} ||\Pi_K(\mathcal{I}(v)) - v||_{L^2(B_K \cap \Omega)}^2 \right), $$

and for any $v \in H^1(\Omega)$

$$ s_{h,0}(\mathcal{I}(v), \mathcal{I}(v)) \le CS_0^* \left( \sum_{K \in \mathcal{T}} ||\mathcal{M}_K(\mathcal{I}(v)) - v||_{L^2(B_K \cap \Omega)}^2 \right), $$

where the constants $C > 0$ only depend on $\Omega$ and the quality parameters of the geometry and the quadrature family. The same results hold replacing $\mathcal{I}$ by $\mathcal{I}^0$ if $v \in H_0^1(\Omega)$. The novelty in the approximation of the advective term comes from the additional stabilization term

$$ s_{h,1}(U, V) = \frac{1}{2} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} |\alpha_\sigma \cdot \eta_{K,\sigma}| (u_\sigma - \mathcal{M}_K(U))(v_\sigma - \mathcal{M}_K(V)) $$

which is easily bounded by:

$$ s_{h,1}(\mathcal{I}(v), \mathcal{I}(v)) \le Ch^2 ||\alpha||_{L^\infty(\Omega)} \theta_\Pi ||\mathcal{I}(v)||_X^2 \quad \text{where} \quad \theta_\Pi = \max_{K \in \mathcal{T}} \max_{\sigma \in \mathcal{F}_K} m_K^{-1} h_K |\eta_{K,\sigma}| $$
where $\theta_{\Pi}$ is one of the quality parameters of the geometry. This gives the consistency of the advective stabilization term as we know from [11] that $||\mathcal{I}(v)||_X$ is controlled by $||v||_{H^1(\Omega)}$.

At this point using the interpolation operator $\mathcal{I}$ one can start to follow the usual lines of error estimation for variational problems, applying the discrete bilinear form to the exact solution tested against any $V \in X_N^0$ and taking the difference with the discrete bilinear form applied to the discrete solution. Doing so one obtains:

$$ 
\begin{aligned}
a_{\Lambda,h}(\mathcal{I}(u) - U, V) + c_{\alpha,h}(\mathcal{I}(u) - U, V) + r_{\beta,h}(\mathcal{I}(u) - U, V) \\
&= \int_{\Omega} \Lambda \nabla u \nabla_{\mathcal{T}}(V) + \frac{1}{2} \int_{\Omega} \boldsymbol{\alpha} \cdot \nabla u \Pi_{\mathcal{T}}(V) - \frac{1}{2} \int_{\Omega} u(\boldsymbol{\alpha}_{\mathcal{F}} \cdot \nabla_{\mathcal{T}})(V) + \frac{1}{2} \int_{\Omega} (\operatorname{div} \boldsymbol{\alpha} + \beta) u \Pi_{\mathcal{T}}(V) \\
&\quad - \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} \zeta_K f \mathcal{M}_K(V) + \mathcal{O}_h \\
&= \int_{\Omega} \Lambda \nabla u \nabla_{\mathcal{T}}(V) + \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} \zeta_K \operatorname{div}(\Lambda \nabla u) \mathcal{M}_K^d(V) \\
&\quad - \frac{1}{2} \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} \zeta_K \operatorname{div}(\boldsymbol{\alpha} u) \mathcal{M}_K(V) - \frac{1}{2} \int_{\Omega} u(\boldsymbol{\alpha}_{\mathcal{F}} \cdot \nabla_{\mathcal{T}})(V) + \mathcal{O}_h,
\end{aligned}
$$

where following [13] $\zeta_{\mathcal{T}} = (\zeta_K)_{K \in \mathcal{T}}$ is a family of functions that satisfies for any $K \in \mathcal{T}$, $\zeta_K \in L^{\infty}(\mathbb{R}^d)$ and:

$$ \int_{\Omega} \zeta_K = m_K, \quad \text{for a.e } \boldsymbol{x} \in \Omega, \quad \operatorname{supp} \zeta_K \subset B_K, \quad M_{\zeta} = \max_{K \in \mathcal{T}} ||\zeta_K||_{L^{\infty}(\Omega)}, $$

and the residual term $\mathcal{O}_h$ is bounded by $\mathcal{O}_h \le C \mathcal{E}_h^c(u, \Lambda, \beta, \boldsymbol{\alpha}, f)$ where the constant $C > 0$ depends on $\Omega$ and the quality parameters but not on $h$ and where the strong consistency error $\mathcal{E}_h^c(u, \Lambda, \beta, \boldsymbol{\alpha}, f)$ is given by:

$$ 
\begin{aligned}
\mathcal{E}_h^c(u, \Lambda, \beta, \boldsymbol{\alpha}, f) = & \left( \sum_{K \in \mathcal{T}} ||\nabla \Pi_K(\mathcal{I}^0(u)) - \nabla u||_{L^2(B_K \cap \Omega)}^2 \right)^{1/2} + \left( \sum_{K \in \mathcal{T}} ||\mathcal{M}_K(\mathcal{I}^0(u)) - u||_{L^2(B_K \cap \Omega)}^2 \right)^{1/2} \\
& + \left( \sum_{K \in \mathcal{T}} h_K^{-2} ||\Pi_K(\mathcal{I}^0(u)) - u||_{L^2(B_K \cap \Omega)}^2 \right)^{1/2} + h^p + \sup_{K \in \mathcal{T}_h} ||\beta - \beta_K||_{L^\infty(B_K \cap \Omega)} + \sup_{K \in \mathcal{T}_h} ||\Lambda - \Lambda_K||_{L^\infty(B_K \cap \Omega)} \\
& + \left( \sum_{K \in \mathcal{T}} ||\mathcal{G}_K(\mathcal{D}_{\mathcal{F}_K}(\boldsymbol{\alpha}), \mathcal{I}_K(u)) - \boldsymbol{\alpha} \cdot \nabla u||_{L^2(B_K \cap \Omega)}^2 \right)^{1/2} + h ||v||_{H^1(\Omega)} + \sup_{K \in \mathcal{T}} ||\mathcal{D} \mathcal{I} \mathcal{V}_K(\boldsymbol{\alpha}_{\mathcal{F}_K}) - \operatorname{div} \boldsymbol{\alpha}||_{L^\infty(B_K)} \\
& + \left( \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} |\zeta_K| |f_K - f|^2 \right)^{1/2}.
\end{aligned}
$$

The strong consistency of the operators $\mathcal{M}_K, \nabla_K$ and $G_K(\boldsymbol{\alpha}_{\mathcal{F}_K}, \cdot)$ with respect to the identity operator, the gradient operator and the advective derivative operator in $B_K$ directly follow from Taylor's expansion and the consistency properties of the network geometry. This then allows under suitable regularity assumptions to bound $\mathcal{E}_h^c(u, \Lambda, \beta, \boldsymbol{\alpha}, f)$ by a constant times $h$. The tricky part for the NEM is precisely to estimate the terms:

$$ W^G(\psi_T, \zeta_T, \Lambda \nabla u, V) = \int_{\Omega} \Lambda \nabla u \cdot \nabla_T(V) + \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} \zeta_K \operatorname{div}(\Lambda \nabla u) \mathcal{M}_K^d(V), $$

and

$$ W^A(\psi_T, \zeta_T, \boldsymbol{\alpha}, u, V) = \frac{1}{2} \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} \zeta_K \operatorname{div}(\boldsymbol{\alpha} u) \mathcal{M}_K(V) + \frac{1}{2} \int_{\Omega} u(\boldsymbol{\alpha}_{\mathcal{F}} \cdot \nabla_T)(V). $$

This corresponds to the additional “limit conformity” or “weak consistency” of the discrete gradient and of the discrete advective derivative that are required to ensure the convergence of discrete variational problems.
Reproducing the proof of [13] for the non conservative case, we easily obtain the following weak consistency estimate on the discrete gradient, where once again the constant $C > 0$ only depend on $\Omega$ and the quality parameters of the geometry and the quadrature family:

$$|W^G(\psi_T, \zeta_T, \Lambda \nabla u, V)| \le C \mathcal{E}_h^G(\zeta_T, \Lambda \nabla u) ||V||_X,$$

where

$$\mathcal{E}_h^G(\zeta_T, \Lambda \nabla u) = \left( \sum_{K \in \mathcal{T}} ||\mathcal{I}_K(\Lambda \nabla u) - \Lambda \nabla u||_{L^2(B_K \cap \Omega)^d}^2 \right)^{1/2} + \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} ||\mathcal{I}_\sigma(\Lambda \nabla u) - \Lambda \nabla u||_{L^2(B_K \cap \Omega)^d}^2 \right)^{1/2} + h^p \\ + \left( \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} |\zeta_K| |\mathcal{D} \mathcal{I} \mathcal{V}_K(\mathcal{I}(\Lambda \nabla u)) - \text{div}(\Lambda \nabla u)|^2 \right)^{1/2},$$

and

$$\mathcal{I}_K(\Lambda \nabla u) = \frac{1}{|B_K \cap \Omega|} \int_{B_K \cap \Omega} \Lambda \nabla u \quad \forall K \in \mathcal{T}.$$

An immediate adaptation leads a similar weak consistency estimate for the advective derivative:

$$|W^A(\psi_T, \zeta_T, \alpha, u, V)| \le C \mathcal{E}_h^A(\zeta_T, \alpha, u) ||V||_X,$$

where

$$\mathcal{E}_h^A(\zeta_T, \alpha, u) = \left( \sum_{K \in \mathcal{T}} ||\mathcal{I}_K(u) - u||_{L^2(B_K \cap \Omega)}^2 \right)^{1/2} + \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{F}_K} ||\mathcal{I}_\sigma(u) - u||_{L^2(B_K \cap \Omega)}^2 \right)^{1/2} + h^p \\ + \left( \sum_{K \in \mathcal{T}} \int_{B_K \cap \Omega} |\zeta_K| |\mathcal{D} \mathcal{I} \mathcal{V}_K(\mathcal{I}_K(\alpha u)) - \text{div} \alpha|^2 \right)^{1/2}.$$

Combining the strong and weak consistency results and using the coercivity and continuity of the global bilinear form, one establishes that:

$$||U - \mathcal{I}^0(u)||_X \le C \mathcal{E}_h(u, \Lambda, \beta, \alpha, f),$$

and

$$||u - \Pi_T(U)||_{L^2(\Omega)} \le C \mathcal{E}_h(u, \Lambda, \beta, \alpha, f) \quad \text{and} \quad ||\nabla u - \nabla_T(U)||_{L^2(\Omega)} \le C \mathcal{E}_h(u, \Lambda, \beta, \alpha, f),$$

where the constant $C > 0$ depends on $\Omega$ and the quality parameters but not on $h$, and where the abstract error is given by:

$$\mathcal{E}_h(u, \Lambda, \beta, \alpha, f) = \mathcal{E}_h^C(u, \Lambda, \beta, \alpha, f) + \mathcal{E}_h^G(\zeta_T, \Lambda \nabla u) + \mathcal{E}_h^A(\zeta_T, \alpha, u).$$

From this result, one finally deduces that for any $K \in \mathcal{T}$:

$$||u - \mathcal{M}_K(U)||_{L^2(B_K \cap \Omega)} \le C \mathcal{E}_h(u, \Lambda, \beta, \alpha, f) \quad \text{and} \quad ||\nabla u - \nabla_K(U)||_{L^2(B_K \cap \Omega)^d} \le C \mathcal{E}_h(u, \Lambda, \beta, \alpha, f).$$

Assuming that $u, \beta, \Lambda$ and $\alpha$ are regular enough, it is now straightforward to bound the total error $\mathcal{E}_h(u, \Lambda, \beta, \alpha, f)$ by a constant times $h$.
REFERENCES

[1] I. Aavatsmark, T. Barkve, O. Boe, and T. Mannseth, Discretization on unstructured grids for inhomogeneous, anisotropic media Part I: derivation of the methods. *SIAM J. Sci. Comput.* **19** (1998) 1700–1716.

[2] I. Aavatsmark, T. Barkve, O. Boe, and T. Mannseth. Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media. *J. Comput. Phys.* **127** (1996) 2–14.

[3] I. Aavatsmark, T. Barkve, O. Boe, and T. Mannseth. Discretization on unstructured grids for inhomogeneous, anisotropic media. Part II: discussion and numerical results. *SIAM J. Sci. Comput.* **19** (1998) 1717–1736.

[4] L. Agélas and R. Masson. Convergence of finite volume mpfa o type schemes for heterogeneous anisotropic diffusion problems on general meshes. *C.R. Acad. Paris, Ser. I* **346** (2008).

[5] L. Agélas, D.A. Di Pietro, and J. Droniou. The g method for heterogeneous anisotropic diffusion on general meshes. *ESAIM: Math. Model. Numer. Anal.* **11** (2010) 597–625.

[6] B. Ahmad, A. Alsaedi, F. Brezzi, L.D. Marini, and A. Russo. Equivalent projectors for virtual element methods. *Comput. Math. Appl.* **66** (2013) 376–391.

[7] I. Babuška, U. Banerjee, and J.E. Osborn. Survey of meshless and generalized finite element methods: a unified approach. *Acta Numer.* **12** (2003) 1–125.

[8] K. Brenner, M. Groza, C. Guichard, and R. Masson. Vertex approximate gradient scheme for hybrid dimensional two-phase darcy flows in fractured porous media. *Finite Volumes for Complex Applications VII-Elliptic, Parabolic and Hyperbolic Problems, Springer Proceedings in Mathematics & Statistics*. Vol. 78 (2014) 507–515.

[9] J.-S. Chen, M. Hillman, and S.-W Chi. Meshfree methods: progress made after 20 years. *J. Eng. Mech.* **143** (2017)

[10] J. Coatléven. A virtual volume method for heterogeneous and anisotropic diffusion-reaction problems on general meshes. *ESAIM: Math. Model. Numeri. Anal.* **51** (2017) 797–824.

[11] J. Coatléven. Basic convergence theory for the network element method. *ESAIM: Math. Model. Numer. Anal.* **55** (2021) 2503–2533.

[12] J. Coatléven. Principles of a network element method. *J. Comput. Phys.* **433** (2021). Art. 110197.

[13] J. Coatléven. A network element method for heterogeneous and anisotropic diffusion-reaction problems. *J. Comput. Phys.* **470** (2022) Art. 111597.

[14] J. Coatléven. On network and geometry generation for the network element method. *preprint*, 2022.

[15] L. Beirao da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L.D. Marini, and A. Russo. Basic principles of virtual element methods. *Math. Models Methods Appl. Sci.* **23** (2013) 199–214.

[16] D. A. Di Pietro and J. Droniou. *The Hybrid High-Order Method for Polytopal Meshes, Design, Analysis, and Applications*. Springer Nature Switzerland (2020).

[17] O. Diyankov. *Uncertain grid method for numerical solution of pdes*. Technical report, NeurOK Software (2008).

[18] J. Droniou. Finite volume schemes for diffusion equations: introduction to and review of modern methods. *M3AS* **24** (2014) 1575–1619. Special edition "P.D.E. Discretizations on Polygonal Meshes".

[19] J. Droniou and R. Eymard. A mixed finite volume scheme for anisotropic diffusion problems on any grid. *Numer. Math.* **105** (2006) 35–71.

[20] P. Durand-Riard, L. Salles, M. Ford, G. Caumon, and J. Pellerin. Understanding the evolution of syn-depositional folds: coupling decompaction and 3d sequential restoration. *Mari. Pet. Geol.* **28** (2011) Art. 1530e1539.

[21] R. Eymard, T. Gallouët, and R. Herbin. A new finite volume scheme for anisotropic diffusion problems on general grids: convergence analysis. *C.R. Math. Acad. Sci. Paris* **344** (2007) 403–406.

[22] R. Eymard, T. Gallouët, and R. Herbin. Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes sushi: a scheme using stabilisation and hybrid interfaces. *IMA J. Numer. Anal.* **30** (2010) 1009–1043.

[23] R. Eymard, C. Guichard, and R. Herbin. Small-stencil 3d schemes for diffusive flows in porous media. *ESAIM: Math. Model. Numer. Anal.* **46** (2011) 265–290.

[24] R. Eymard, C. Guichard, R. Herbin, and R. Masson. Vertex-centred discretization of multiphase compositional darcy flows on general meshes. *Comput. Geosci.* **16** (2012) 987–1005.

[25] R. Eymard, C. Guichard, R. Herbin, and R. Masson. Vertex centred discretization of two-phase darcy flows on general meshes. *ESAIM: Proc.* **35** (2012) 59–78.

[26] A. Katz and A. Jameson. A meshless volume scheme. In: *Proc. of 19th AIAA Computational Fluid Dynamics, Fluid Dynamics and Co-located Conferences*. American Institute of Aeronautics and Astronautics (2009) 2009–3534.

[27] A. I. Kauerauf and T. Hantschel. *Fundamentals of Basin and Petroleum Systems Modeling*. Springer Berlin, Heidelberg (2009).

[28] J.M. Melenk. *On Approximation in Meshless Methods*. Springer Berlin Heidelberg, Berlin, Heidelberg (2005) 65–141.

[29] L. Ran, H. Borouchaki, A. Benali, and C. Bennis. Hex-dominant mesh generation for basin modeling with complex geometry. *IOP Conf. Ser.: Mater. Sci. Eng.* **10** (2010) 012085.

[30] V. Shankar, G. B. Wright, and A. L. Fogelson. An efficient high-order meshless method for advection-diffusion equations on time-varying irregular domains. *J. Comput. Phys.* **445** (2021) Art. 110633.

[31] J. Slak and G. Kosec. On generation of node distributions for meshless pde discretizations. *SIAM J. Sci. Comput.* **41** (2019) A3202–A3229.

[32] E.M. Stein. *Singular integrals and differentiability properties of functions*. Princeton University Press (1970).

[33] A. Tertois and J. Mallet. Restoration of complex three-dimensional structural models based on the mathematical geochron framework. In: *81st EAGE Conference and Exhibition 2019*. Vol. 2019 (2019) 1–5.
[34] N. Trask, P. Bochev, and M. Perego. A conservative, consistent, and scalable mesh-free mimetic method. *J. Comput. Phys.* **409** (2020) 109–187.

[35] N. Trask, M. Perego, and P. Bochev. A high-order staggered meshless method for elliptic problems. *SIAM J. Sci. Comput.* **39** (2017) 479–502.

[36] F. Wellmann and G. Caumon. Chapter one – 3-d structural geological models: Concepts, methods, and uncertainties. In: Vol. 59 of *Advances in Geophysics*. Elsevier (2018) 1–121.

[37] E. Kwan yu Chiu, Q. Wang, R. Hu, and A. Jameson. A conservative mesh-free scheme and generalized framework for conservation laws. *SIAM J. Sci. Comput.* **34** (2012) 2896–2916.

<figure><img src="image_15.png" /></figure>

Please help to maintain this journal in open access!

This journal is currently published in open access under the **Subscribe to Open model (S2O)**. We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.