# Finite Volumes for the Stefan-Maxwell Cross-Diffusion System

Clément Cancès, Virginie Ehrlacher, Laurent Monasse

► To cite this version:

Clément Cancès, Virginie Ehrlacher, Laurent Monasse. Finite Volumes for the Stefan-Maxwell Cross-Diffusion System. IMA Journal of Numerical Analysis, In press, 10.1093/imanum/drad032. hal-02902672v2

HAL Id: hal-02902672

https://hal.science/hal-02902672v2

Submitted on 9 Jan 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
## FINITE VOLUMES FOR THE STEFAN-MAXWELL CROSS-DIFFUSION SYSTEM

CLÉMENT CANCÈS, VIRGINIE EHRLACHER, AND LAURENT MONASSE

**ABSTRACT.** The aim of this work is to propose a provably convergent finite volume scheme for the so-called Stefan-Maxwell model, which describes the evolution of the composition of a multi-component mixture and reads as a cross-diffusion system. The scheme proposed here relies on a two-point flux approximation, and preserves at the discrete level some fundamental theoretical properties of the continuous models, namely the non-negativity of the solutions, the conservation of mass and the preservation of the volume-filling constraints. In addition, the scheme satisfies a discrete entropy-entropy dissipation relation, very close to the relation which holds at the continuous level. In this article, we present this scheme together with its numerical analysis, and finally illustrate its behaviour with some numerical results.

### 1. THE STEFAN-MAXWELL MODEL

The aim of this section is to present the so-called Stefan-Maxwell model, which is introduced in Section 1.1. Its key mathematical properties are summarized in Section 1.2. In particular, an entropy-entropy dissipation inequality holds for this system and is formally derived in Section 1.3.

1.1. **Presentation of the model.** The Maxwell-Stefan equations describe the evolution of the composition of a multicomponent mixture via diffusive transport [45, 49]. This model is used in various applications like sedimentation, dialysis, electrolysis, ion exchange, ultrafiltration, and respiratory airways [54].

We are interested in the evolution of the composition of a mixture of $n \in \mathbb{N}^*$ species, which is described by the volume fractions $u = (u_1, \dots, u_n)$, where $u_i$ denotes the volume fraction of the $i^{th}$ species for all $1 \le i \le n$. The spatial domain occupied by the mixture is represented by an open, connected, bounded, and polyhedral subset $\Omega$ of $\mathbb{R}^d$. Let $T > 0$ denote some arbitrary final time.

For all $1 \le i \ne j \le n$, let $c_{ij} = c_{ji} > 0$ be some positive real numbers. The coefficient $c_{ij}$ can be interpreted as the inverse of the inter-species diffusion coefficient between the $i^{th}$ and $j^{th}$ species. For all $v := (v_1, \dots, v_n) \in \mathbb{R}_+^n$, we denote by $A(v) := (A_{ij}(v))_{1 \le i,j \le n}$ the matrix defined by

$$ (1) \qquad A_{ii}(v) := \sum_{1 \le j \ne i \le n} c_{ij}v_j, \qquad A_{ij}(v) := -c_{ij}v_i. $$

In the Stefan-Maxwell model, the evolution of the composition of the mixture is prescribed by the following system of partial differential equations:

$$ (2) \qquad \partial_t u_i + \text{div} J_i = 0, \quad \forall 1 \le i \le n, $$
where the set of fluxes $J := (J_i)_{1\le i\le n}$ is solution to the set of equations

$$ (3) \qquad \nabla u_i + \sum_{j=1}^{n} A_{ij}(u) J_j = 0, \quad \forall 1 \le i \le n, $$

$$ (4) \qquad \sum_{i=1}^{n} J_i = 0. $$

For any vectors $v := (v_i)_{1\le i\le n}$, $w := (w_i)_{1\le i\le n} \in \mathbb{R}^n$, we denote by $\langle v, w \rangle := \sum_{i=1}^n v_i w_i$ the canonical scalar product of $v, w$ in $\mathbb{R}^n$, while the canonical scalar product of vectors $F, G \in \mathbb{R}^d$ is denoted by $F \cdot G$. Equations (3) and (4) can then be rewritten in the more compact form

$$ (5) \qquad \nabla u + A(u)J = 0, $$

$$ (6) \qquad \langle \mathbf{1}, J \rangle = 0, $$

where $\mathbf{1} := (1, 1, \dots, 1) \in \mathbb{R}^n$. We refer the reader to Appendix A of [41] and [12] for the derivation of the model (2)-(3)-(4).

The system is complemented with no-flux boundary conditions

$$ (7) \qquad J_i \cdot \nu = 0 \text{ on } \partial\Omega, \quad \text{for all } 1 \le i \le n, $$

and a measurable initial condition $u^0 = (u_1^0, \dots, u_n^0)$ which satisfies

$$ (8) \qquad \forall 1 \le i \le n, \quad u_i^0 \ge 0 \quad \text{and} \quad \sum_{i=1}^{n} u_i^0 = 1 \quad \text{in } \Omega. $$

In other words, denoting by

$$ \mathcal{A} := \{ v \in \mathbb{R}_+^n, \langle \mathbf{1}, v \rangle = 1 \}, $$

we assume that $u^0 \in L^\infty(\Omega; \mathcal{A})$. Let us also assume in addition that

$$ (9) \qquad \forall 1 \le i \le n, \quad M_i := \int_{\Omega} u_i^0 > 0, $$

i.e. that each of the different species is initially present in the mixture. We denote by $M = (M_i)_{1\le i\le n} \in (\mathbb{R}_+^*)$ the vector of masses. Since $u_0 \in L^\infty(\Omega; \mathcal{A})$, one has $\langle \mathbf{1}, M \rangle = m_\Omega$ where $m_\Omega$ stands for the Lebesgue measure of $\Omega$.

The mathematical analysis of the Stefan-Maxwell model is quite recent [33, 10, 11, 41]. The first existence result of global weak solutions to the Stefan-Maxwell problem for general initial data and number of chemical species was proved in [41].

Motivated by the results of [41], we introduce here the notion of weak solution to the Stefan-Maxwell system of equations, which is used in our analysis. In what follows, we denote by $Q_T = (0, T) \times \Omega$, and by

$$ \mathcal{V}_\lambda = \left\{ v = (v_1, \dots, v_n) \left| \sum_i v_i = \langle \mathbf{1}, v \rangle = \lambda \right. \right\}, \quad \lambda \in \mathbb{R}. $$

In particular, $\mathcal{A} = \mathcal{V}_1 \cap (\mathbb{R}_+)^n$, $J \in (\mathcal{V}_0)^d$, and $M \in \mathcal{V}_{m_\Omega}$.

**Definition 1.1.** A weak solution $(u, J)$ to (2)-(5)-(6) corresponding to the initial profile $u^0 \in L^\infty(\Omega; \mathcal{A})$ is a pair $(u, J)$ such that $u \in L^\infty(Q_T; \mathcal{A}) \cap L^2((0, T); H^1(\Omega))^n$
and $\nabla\sqrt{u} \in L^2(Q_T)^{n \times d}$, such that $J \in L^2(Q_T; (V_0)^d)$ satisfies (5), and such that, for all $\phi \in C_c^\infty([0, T) \times \bar{\Omega})^n$,

$$ (10) \qquad \iint_{Q_T} \langle u, \partial_t \phi \rangle + \int_{\Omega} \langle u^0, \phi(0, \cdot) \rangle + \iint_{Q_T} \sum_{i=1}^{n} J_i \cdot \nabla \phi_i = 0. $$

To the best of our knowledge, the uniqueness of weak solutions are still an open problem. However, such weak solutions are unique provided one strong solution exists due to some recent weak-strong uniqueness result [37]. Let us finally mention the partial regularity result [14] which shows that the weak solutions are regular up to a small set.

## 1.2. Key mathematical properties of the model.
In this section, we exhibit some key mathematical properties of the model, which were proved in [41], and that we wish to preserve at the discrete level in the numerical scheme.

First, the total mass of each species is conserved, i.e, for all $1 \le i \le n$ and $t > 0$,

$$ (11) \qquad \int_{\Omega} u_i(t, x) \, dx = \int_{\Omega} u_i^0(x) \, dx. $$

This follows directly from the local conservation property (2) and the no-flux boundary conditions across $\partial\Omega$.

Second, the volume fractions remain non-negative, i.e.,

$$ (12) \qquad \forall 1 \le i \le n, \quad u_i(t, x) \ge 0, \text{ for almost all } (t, x) \in Q_T. $$

Third, the condition (4) together with (2) implies that $\partial_t \langle \mathbf{1}, u \rangle = 0$. Therefore, the volume filling property (8) satisfied by the initial condition is satisfied for all $t > 0$:

$$ (13) \qquad \sum_{i=1}^{n} u_i(t, x) = 1 \quad \text{for almost all } (t, x) \in Q_T. $$

In particular, $u \in L^\infty(Q_T; \mathcal{A})$.

Lastly, an entropy-entropy dissipation relation, which is formally derived in Section 1.3, holds for this system, so that the functional

$$ E: \begin{cases} L^\infty(\Omega, \mathcal{A}) & \to \mathbb{R} \\ u := (u_1, \dots, u_n) & \mapsto \int_{\Omega} \sum_{i=1}^{n} u_i \log u_i \end{cases} $$

is a Lyapunov function for the Stefan-Maxwell system. More precisely, it holds that

$$ (14) \qquad \frac{d}{dt} E(u(t)) + \frac{\alpha}{2} \int_{\Omega} \sum_{i=1}^{n} |\nabla \sqrt{u_i}|^2 + \frac{c^*}{2} \int_{\Omega} \sum_{i=1}^{n} |J_i|^2 \le 0, $$

for some positive constants $\alpha, c^* > 0$ whose definitions are made precise in the next section.

## 1.3. Continuous entropy estimate.
We formally derive here the entropy-entropy dissipation inequality (14) which holds for the continuous system and was rigorously proved in [41]. For the formal calculations to hold, we make the simplifying assumption in this Section that the solution $u$ to the Stefan-Maxwell model satisfies

$$ (15) \qquad \forall 1 \le i \le n, \quad u_i(t, x) > 0 \quad \text{and} \quad \sum_{i=1}^{n} u_i(t, x) = 1 \quad \text{a.e. in } Q_T, $$

and that the solution enjoys enough regularity to justify the calculations.
To present the entropy-entropy dissipation inequality which holds for the Stefan-Maxwell model, we need to introduce some additional notation. Denote by

$$c^* = \min_{1 \le i \ne j \le n} c_{ij} > 0,$$

then for all $1 \le i \ne j \le n$, we define

$$\bar{c}_{ij} := c_{ij} - c^* \quad \text{and} \quad \bar{c} := \max_{1 \le i \ne j \le n} \bar{c}_{ij}.$$

Let us point out that $\bar{c}_{ij} \ge 0$ for all $1 \le i \ne j \le n$ (and thus $\bar{c} \ge 0$).

Let I denote the $n \times n$ identity matrix. For all $v \in \mathbb{R}^n$, we introduce $\bar{A}(v) := (\bar{A}_{ij}(v))_{1 \le i,j \le n}$ and $C(v) := (C_{ij}(v))_{1 \le i,j \le n}$ the matrices respectively defined as follows: for all $1 \le i, j \le n$,

$$(16) \qquad \bar{A}_{ii}(v) := \sum_{1 \le j \ne i \le n} \bar{c}_{ij}v_j, \quad \bar{A}_{ij}(v) := -\bar{c}_{ij}v_i \quad \text{and} \quad C_{ij}(v) := v_i.$$

It then holds that for all $v := (v_1, \dots, v_n) \in (\mathbb{R}_+)^n$,

$$(17) \qquad A(v) = c^* \langle 1, v \rangle I - c^* C(v) + \bar{A}(v).$$

In particular, if $u \in \mathbb{R}_+^n$ satisfies $\langle 1, u \rangle = 1$, then

$$(18) \qquad A(u) = c^* I - c^* C(u) + \bar{A}(u).$$

One easily deduces from the particular form (16) of the matrix $\bar{A}(v)$ that

$$(19) \qquad \mathrm{Span}\{v\} \subset \mathrm{Ker}(\bar{A}(v)), \quad \mathrm{Ran}(\bar{A}(v)) \subset \mathcal{V}_0, \quad \forall v \in \mathbb{R}^n.$$

It has been established in [41] that equalities instead of mere inclusions hold in (19) if one replaces $\bar{A}(v)$ by $A(v)$ and one considers $v$ with positive components, i.e.,

$$(20) \qquad \mathrm{Span}\{v\} = \mathrm{Ker}(A(v)), \quad \mathrm{Ran}(A(v)) = \mathcal{V}_0, \quad \forall v \in (\mathbb{R}_+^*)^n.$$

This property is intensively used in the convergence study of [41]. Provided (15) holds, (20) shows that there exists a unique solution $J(t, x)$ to (5)-(6) for almost all $(t, x) \in (0, T) \times \Omega$, since $\nabla u \in (\mathcal{V}_0)^d$. Besides, using (18), it holds that $J$ is a solution to (5)-(6) if and only if it is the unique solution to

$$(21) \qquad \nabla u + c^* J + \bar{A}(u) J = 0, \quad \forall 1 \le i \le n,$$

$$(22) \qquad \langle 1, J \rangle = 0,$$

since $\langle 1, u \rangle = 1$ and since the condition $\langle 1, J \rangle = 0$ implies that $C(u)J = 0$.

For all $v := (v_1, \dots, v_n) \in (\mathbb{R}_+^*)^n$, we denote by $M(v) := \mathrm{diag}(v_1, \dots, v_n)$ the $n \times n$ diagonal matrix whose $i^{th}$ diagonal entry is given by $v_i$ for all $1 \le i \le n$. Then, the following lemma, which is central in our analysis, holds.

**Lemma 1.2.** Let $v := (v_1, \dots, v_n) \in (\mathbb{R}_+^*)^n$, such that $\langle 1, v \rangle \le 1$. Then, it holds that $\bar{B}(v) := M^{-1}(v)\bar{A}(v)$ is a symmetric semi-definite non-negative matrix such that

$$(23) \qquad M^{-1}(v)\bar{A}(v) \le 2\bar{c}M^{-1}(v),$$

in the sense of symmetric matrices.
Proof. Let $v := (v_1, \cdots, v_n) \in (\mathbb{R}_+^*)^n$ and $\bar{B}(v) := M^{-1}(v)\bar{A}(v)$. Denoting by $(\bar{B}_{ij}(v))_{1\le i,j\le n}$ the different components of $\bar{B}(v)$, a direct calculation shows that for all $1 \le i, j \le n$,

$$ \bar{B}_{ij}(v) := -\bar{c}_{ij} \text{ if } i \neq j \quad \text{and} \quad \bar{B}_{ii}(v) = \sum_{1 \le j \neq i \le n} \bar{c}_{ij} \frac{v_j}{v_i}, $$

hence the symmetry of the matrix $\bar{B}(v)$. Let $\xi := (\xi_i)_{1\le i\le n} \in \mathbb{R}^n$. Using the fact that $\bar{c}_{ij} = \bar{c}_{ji}$ for all $1 \le i \neq j \le n$, it holds that

$$ \begin{aligned} \xi^T \bar{B}(v) \xi &= \sum_{1 \le j \neq i \le n} \bar{c}_{ij} \left( \frac{v_j}{v_i} \xi_i^2 - \xi_i \xi_j \right) \\ &= \frac{1}{2} \sum_{1 \le j \neq i \le n} \bar{c}_{ij} \left( \frac{v_j}{v_i} \xi_i^2 + \frac{v_i}{v_j} \xi_j^2 - 2\xi_i \xi_j \right) \\ &= \frac{1}{2} \sum_{1 \le j \neq i \le n} \bar{c}_{ij} \left( \sqrt{\frac{v_j}{v_i}} \xi_i - \sqrt{\frac{v_i}{v_j}} \xi_j \right)^2 \ge 0. \end{aligned} $$

Hence the non-negativity of the matrix $\bar{B}(v)$. Using now the elementary inequality $(a-b)^2 \le 2a^2 + 2b^2$ together with the fact that $\langle 1, v \rangle \le 1$ for all $1 \le i \le n$, we obtain that

$$ \begin{aligned} \xi^T \bar{B}(v) \xi &= \frac{1}{2} \sum_{1 \le j \neq i \le n} \bar{c}_{ij} \left( \sqrt{\frac{v_j}{v_i}} \xi_i - \sqrt{\frac{v_i}{v_j}} \xi_j \right)^2 \\ &\le \sum_{1 \le j \neq i \le n} \bar{c}_{ij} \left( \frac{v_j}{v_i} \xi_i^2 + \frac{v_i}{v_j} \xi_j^2 \right) \\ &\le \bar{c} \sum_{1 \le j \neq i \le n} \left( \frac{1}{v_i} \xi_i^2 + \frac{1}{v_j} \xi_j^2 \right) \le 2\bar{c} \xi^T M^{-1}(v) \xi. \end{aligned} $$

Hence the desired result.

□

We are now in position to write the (formal) entropy-entropy dissipation inequality which holds on the continuous level for the Stefan-Maxwell model. For all $1 \le i \le n$, let $w_i = \log u_i$ and $w = (w_i)_{1\le i\le n}$. Then, it holds that $\nabla u = M(u)\nabla w$ which implies that

$$ (24) \qquad \nabla w = -M^{-1}(u)A(u)J = -(c^*M^{-1}(u) + M^{-1}(u)\bar{A}(u))J. $$

Since $M^{-1}(u)$ is symmetric definite positive while $M^{-1}(u)\bar{A}(u)$ is symmetric non-negative, it holds that $c^*M^{-1}(u) + M^{-1}(u)\bar{A}(u)$ is an invertible matrix so that $J = -(c^*M^{-1}(u) + M^{-1}(u)\bar{A}(u))^{-1}\nabla w$. This yields that

$$ (25) \qquad \frac{d}{dt} E(u(t)) = \int_{\Omega} \sum_{i=1}^{n} \partial_t u_i w_i \stackrel{(2)}{=} - \int_{\Omega} \sum_{i=1}^{n} \mathrm{div} J_i w_i \stackrel{(7)}{=} \int_{\Omega} J \cdot \nabla w. $$

Using (24), the last term in the above equality can be rewritten of two different manners:

$$ (26) \qquad \int_{\Omega} J \cdot \nabla w = - \int_{\Omega} J \cdot (c^* M^{-1}(u) + M^{-1}(u) \bar{A}(u)) J $$

$$ (27) \qquad = - \int_{\Omega} \nabla w \cdot (c^* M^{-1}(u) + M^{-1}(u) \bar{A}(u))^{-1} \nabla w. $$
Define the matrix

$$ (28) \quad B(v) := (c^* M^{-1}(v) + M^{-1}(v) \bar{A}(v)), \quad \forall v = (v_i)_{1\le i\le n} \in (\mathbb{R}_+^*)^n. $$

It follows from Lemma 1.2 that the two inequalities

$$ (29) \quad B(v) \ge c^* M^{-1}(v) \ge c^* I, \quad B(v)^{-1} \ge \frac{1}{c^* + 2\bar{c}} M(v), \quad \forall v \in (0, 1]^n, $$

hold in the sense of symmetric matrices. Therefore, we obtain from (26)-(27) that

$$ (30) \quad \int_{\Omega} J \cdot \nabla w \ge -\frac{1}{2(c^* + 2\bar{c})} \int_{\Omega} \nabla w \cdot M(u) \nabla w - \frac{c^*}{2} \int_{\Omega} |J|^2. $$

The first term of the right-hand side can be rewritten by noticing that

$$ \nabla w \cdot M(u) \nabla w = \sum_{i=1}^{n} u_i \nabla \log(u_i) \cdot \nabla \log(u_i) = 4 \sum_{i=1}^{n} |\nabla \sqrt{u_i}|^2. $$

As a consequence, we finally deduce from (25) and (30) that

$$ \frac{d}{dt} E(u(t)) \le -\frac{1}{2} \alpha \int_{\Omega} \sum_{i=1}^{n} |\nabla \sqrt{u_i}|^2 - \frac{1}{2} c^* \int_{\Omega} |J|^2, $$

with

$$ (31) \qquad \alpha := \frac{4}{c^* + 2\bar{c}} > 0. $$

This entropy-entropy dissipation inequality is similar to (14).

**Remark 1.3.** Since the entropy $E$ is bounded on $L^\infty(Q_T; A)$ — it takes its values in $[-m_\Omega \log(n), 0]$ — integrating (14) over $t \in (0, T)$ yields

$$ \iint_{Q_T} |\nabla \sqrt{u}|^2 + \iint_{Q_T} |J|^2 \le C. $$

Moreover, since $u$ is uniformly bounded between 0 and 1, one has

$$ \iint_{Q_T} |\nabla \sqrt{u}|^2 \ge \frac{1}{4} \iint_{Q_T} |\nabla u|^2, $$

so that one gets a control over the $L^2(0, T; H^1(\Omega))$ norm of $u$ and on the $L^2(Q_T)$ norm of $J$. This motivates the weak formulation used in Definition 1.1.

**1.4. Contributions and positionning of the paper.** The goal of this paper is to build and analyze a numerical scheme preserving the properties discussed in the previous section, namely:

* the non-negativity of the concentrations;

* the conservation of mass;

* the preservation of the volume filling constraint;

* the entropy-entropy dissipation relation (14).

The scheme proposed here relies on two-point flux approximation (TPFA) finite volumes [28, 27] and builds on similar ideas as the one introduced in [18] for another family of cross-diffusion systems.

TPFA finite volumes is popular to approximate conservation laws. Unsurprisingly, schemes entering this family of methods have been proposed for the Stefan-Maxwell diffusion problem in [50, 11, 47]. Those schemes yield satisfactory numerical outputs but there is no theoretical guarantee of their convergence so far.
Besides, a finite element scheme is proposed and analysed in [40] for the more complex case where the chemical species are ions inducing a self-consistent electrical potential. The analysis carried out in [40] relies on the very strong assumption that integrals of non-polynomial functions can be computed exactly.

Convergence proofs for finite volume approximations of cross-diffusion systems have been proposed in [3, 1, 21, 17, 18, 42, 23, 46, 32, 31]. Most of the above contributions rely on the entropy-stability of the schemes, which is exploited thanks to the so-called discrete entropy method [22]. This approach is a transposition to the discrete setting of the boundedness-by-entropy method exposed in [38, 39]. The design of entropy stable numerical schemes for diffusion type equations has received an important attention in the last years. Let us mention the contributions [7, 9, 8, 19, 20, 15, 43, 51, 2, 52, 48, 16], this list being non-exhaustive. We mention in particular the recent work [36] where the authors propose an energy stable and positivity-preserving scheme for the Maxwell-Stefan diffusion system, but a convergence analysis of the scheme proposed therein seems out of reach because of the explicit in time treatment of the mobility matrix B.

Let us also mention that finite element methods are also used for the simulation of cross-diffusion systems. We refer the reader to [29, 6, 34] for more details. We would like to highlight the work [13] where the authors propose a space-time Galerkin method which preserves the entropy structure of cross-diffusion systems, including the Stefan-Maxwell system under consideration, as well as the contribution [53] where steady states are computed thanks to finite elements. The coupling with compressible Stokes equation is addressed in [5].

The scheme is presented in Section 2. Our main results are gathered in Section 2.3. Preliminary estimates and existence of a solution to the discretized scheme are proved in Section 3. Convergence of the discretized solution to a weak solution of the continuous model is proved in Section 4. Finally, numerical tests illustrating the behaviour of the method are presented in Section 5.

## 2. THE FINITE-VOLUME SCHEME

### 2.1. Discretization of $(0, T) \times \Omega$.
As already mentioned, our scheme relies on TPFA finite volumes. As explained in [25, 27, 30], this approach appears to be very efficient as soon as the continuous problem to be solved numerically is isotropic and one has the freedom to choose a suitable mesh fulfilling the so-called orthogonality condition [35, 28]. We recall here the definition of such a mesh.

**Definition 2.1.** An admissible mesh of $\Omega$ is a triplet $(\mathcal{T}, \mathcal{E}, (x_K)_{K \in \mathcal{T}})$ such that the following conditions are fulfilled.

(i) Each control volume (or cell) $K \in \mathcal{T}$ is non-empty, open, polyhedral and convex. We assume that

$$K \cap L = \emptyset \text{ if } K, L \in \mathcal{T} \text{ with } K \neq L, \quad \text{while } \bigcup_{K \in \mathcal{T}} \bar{K} = \bar{\Omega}.$$

(ii) Each face $\sigma \in \mathcal{E}$ is closed and is contained in a hyperplane of $\mathbb{R}^d$, with positive $(d-1)$-dimensional Hausdorff (or Lebesgue) measure denoted by $m_\sigma = \mathcal{H}^{d-1}(\sigma) > 0$. We assume that $\mathcal{H}^{d-1}(\sigma \cap \sigma') = 0$ for $\sigma, \sigma' \in \mathcal{E}$ unless $\sigma = \sigma'$. For all $K \in \mathcal{T}$, we assume that there exists a subset $\mathcal{E}_K$ of $\mathcal{E}$ such that $\partial K = \bigcup_{\sigma \in \mathcal{E}_K} \sigma$. Moreover, we suppose that $\bigcup_{K \in \mathcal{T}} \mathcal{E}_K = \mathcal{E}$. Given two distinct control volumes $K, L \in \mathcal{T}$, the intersection $\bar{K} \cap \bar{L}$ either reduces to
a single face $\sigma \in \mathcal{E}$ denoted by $K|L$, or its $(d-1)$-dimensional Hausdorff measure is 0.

(iii) The cell-centers $(x_K)_{K \in \mathcal{T}}$ satisfy $x_K \in K$, and are such that, if $K, L \in \mathcal{T}$ share a face $K|L$, then the vector $x_L - x_K$ is orthogonal to $K|L$.

We denote by $m_K$ the $d$-dimensional Lebesgue measure of the control volume $K$. The set of the faces is partitioned into two subsets: the set $\mathcal{E}_{\text{int}}$ of the interior faces defined by

$$ \mathcal{E}_{\text{int}} = \{ \sigma \in \mathcal{E} \mid \sigma = K|L \text{ for some } K, L \in \mathcal{T} \}, $$

and the set $\mathcal{E}_{\text{ext}} = \mathcal{E} \setminus \mathcal{E}_{\text{int}}$ of the exterior faces defined by $\mathcal{E}_{\text{ext}} = \{\sigma \in \mathcal{E} \mid \sigma \subset \partial\Omega\}$. For a given control volume $K \in \mathcal{T}$, we also define $\mathcal{E}_{K,\text{int}} = \mathcal{E}_K \cap \mathcal{E}_{\text{int}}$ (respectively $\mathcal{E}_{K,\text{ext}} = \mathcal{E}_K \cap \mathcal{E}_{\text{ext}}$) the set of its faces that belong to $\mathcal{E}_{\text{int}}$ (respectively $\mathcal{E}_{\text{ext}}$). For such a face $\sigma \in \mathcal{E}_{K,\text{int}}$, we may write $\sigma = K|L$, meaning that $\sigma = \bar{K} \cap \bar{L}$, where $L \in \mathcal{T}$.

Given $\sigma = K|L \in \mathcal{E}_{\text{int}}$, we let $d_{\sigma} = |x_K - x_L|$ and $\tau_{\sigma} = \frac{m_{\sigma}}{d_{\sigma}}$. We also define $d_{K\sigma} = \text{dist}(x_K, \sigma)$ and $\tau_{K\sigma} = \frac{m_{\sigma}}{d_{K\sigma}}$. Moreover, for all $K \in \mathcal{T}$ and all $\sigma = K|L \in \mathcal{E}_{K,\text{int}}$, we denote by $\nu_{K\sigma} = \frac{x_L - x_K}{d_{K\sigma}}$, the unitary normal to $\sigma$ outward with respect to $K$. The half-diamond cell $\Delta_{K\sigma}$ associated to $K$ and $\sigma$ is defined as the convex hull of $x_K$ and $\sigma$, and we define the diamond cells $\Delta_{\sigma} = \Delta_{K\sigma} \cup \Delta_{L\sigma}$ for $\sigma = K|L$. Then it follows from the an elementary geometrical property that the $(d$-dimensional) Lebesgue measures of $\Delta_{\sigma}$ (resp. $\Delta_{K\sigma}$) are given by

$$ (32) \qquad m_{\Delta_{\sigma}} = \frac{m_{\sigma} d_{\sigma}}{d}, \quad m_{\Delta_{K\sigma}} = \frac{m_{\sigma} d_{K\sigma}}{d}, $$

We finally introduce the size $h_{\mathcal{T}}$ and the regularity $\zeta_{\mathcal{T}}$ (which is assumed to be positive) of a discretization $(\mathcal{T}, \mathcal{E}, (x_K)_{K \in \mathcal{T}})$ of $\Omega$ by setting

$$ h_{\mathcal{T}} = \max_{K \in \mathcal{T}} \text{diam}(K) \quad \text{and} \quad \zeta_{\mathcal{T}} = \min_{K \in \mathcal{T}} \min_{\sigma \in \mathcal{E}_{K,\text{int}}} \frac{d(x_K, \sigma)}{d_{\sigma}}. $$

Concerning the time discretization of $(0, T)$, we consider $P_T \in \mathbb{N}^*$ and an increasing infinite family of times $0 = t_0 < t_1 < \dots < t_{P_T} = T$. We denote by $\Delta t_p = t_p - t_{p-1}$ for $p \in \{1, \dots, P_T\}$, by $\Delta t = (\Delta t_p)_{1 \le p \le P_T}$, and by $h_T = \max_{1 \le p \le P_T} \Delta t_p$. In what follows, we will use boldface notation for mesh-indexed families, typically for elements of $\mathbb{R}^T$, $\mathbb{R}^{\mathcal{E}}$, $(\mathbb{R}^T)^n$, $(\mathbb{R}^{\mathcal{E}})^n$, $(\mathbb{R}^T)^{P_T}$, $(\mathbb{R}^{\mathcal{E}})^{P_T}$ or even $(\mathbb{R}^T)^{n \times P_T}$ and $(\mathbb{R}^{\mathcal{E}})^{n \times P_T}$. One naturally defines discrete $L^2$ scalar products on $\mathbb{R}^T$ and $\mathbb{R}^{d \times \mathcal{E}}$ by setting

$$ \langle \mathbf{u}, \mathbf{v} \rangle_{\mathcal{T}} = \sum_{K \in \mathcal{T}} m_K u_K v_K, \qquad \mathbf{u} = (u_K)_{K \in \mathcal{T}}, \mathbf{v} = (v_K)_{K \in \mathcal{T}} \in \mathbb{R}^T $$

and

$$ \langle \mathbf{F}, \mathbf{G} \rangle_{\mathcal{E}} = \sum_{\sigma \in \mathcal{E}_{\text{int}}} m_{\Delta_{\sigma}} F_{K\sigma} \cdot G_{K\sigma}, \qquad \mathbf{F} = (F_{K\sigma})_{\sigma \in \mathcal{E}}, \mathbf{G} = (G_{K\sigma})_{\sigma \in \mathcal{E}} \in \mathbb{R}^{d \times \mathcal{E}}. $$

Remark that the above definition of $\langle \cdot, \cdot \rangle_{\mathcal{E}}$ does not depend on the choice of the element $K \in \mathcal{T}$ such that $\sigma \in \mathcal{E}_{K,\text{int}}$ for internal edges.
2.2. **Numerical scheme.** The initial data $u^0 \in L^\infty(\Omega; \mathcal{A})$ is discretized into

$$ \boldsymbol{u}^0 = (\boldsymbol{u}_i^0)_{1\le i \le n} \in (\mathbb{R}^T)^n = (u_{i,K}^0)_{K \in \mathcal{T}, 1 \le i \le n}, $$

by setting

$$ (33) \qquad u_{i,K}^0 = \frac{1}{m_K} \int_K u_i^0(x) dx, \quad \forall K \in \mathcal{T}, 1 \le i \le n. $$

Assume that $\boldsymbol{u}^{p-1} = (u_{i,K}^{p-1})_{K \in \mathcal{T}, 1 \le i \le n}$ is given for some $p \ge 1$, then we have to define how to compute the discrete volume fractions $\boldsymbol{u}^p = (u_{i,K}^p)_{K \in \mathcal{T}, 1 \le i \le n}$ and the discrete fluxes $\boldsymbol{J}^p = (J_{i,K\sigma}^p)_{\sigma \in \mathcal{E}, 1 \le i \le n}$.

First, we introduce some notation. Given any discrete scalar field $\boldsymbol{v} = (v_K)_{K \in \mathcal{T}} \in \mathbb{R}^T$, we define for all cell $K \in \mathcal{T}$ and interface $\sigma \in \mathcal{E}_K$ the mirror value $v_{K\sigma}$ of $v_K$ across $\sigma$ by setting:

$$ v_{K\sigma} = \begin{cases} v_L & \text{if } \sigma = K |L \in \mathcal{E}_{\text{int}}, \\ v_K & \text{if } \sigma \in \mathcal{E}_{\text{ext}}. \end{cases} $$

We also define the oriented and absolute jumps of $\boldsymbol{v}$ across any edge by

$$ D_{K\sigma}\boldsymbol{v} = v_{K\sigma} - v_K, \quad \text{and} \quad D_{\sigma}\boldsymbol{v} = |D_{K\sigma}\boldsymbol{v}|, \quad \forall K \in \mathcal{T}, \forall \sigma \in \mathcal{E}_K. $$

Note that in the above definition, for all $\sigma \in \mathcal{E}$, the definition of $D_{\sigma}\boldsymbol{v}$ does not depend on the choice of the element $K \in \mathcal{T}$ such that $\sigma \in \mathcal{E}_K$.

For all $1 \le i \le n$, we also introduce some edge values $u_{i,\sigma}^p$ of the volume fraction $u_i$ for all $\sigma \in \mathcal{E}$. For any $K \in \mathcal{T}$ such that $\sigma \in \mathcal{E}_K$, the definition of $u_{i,\sigma}^p$ makes use of the values $u_{i,K}^p$ and $u_{i,K\sigma}^p$ but is independent of the choice of $K$. As in [18], the edge volume fraction $u_{i,\sigma}^p$ is defined through a logarithmic mean as follows

$$ (34a) \qquad u_{i,\sigma}^p = \begin{cases} 0 & \text{if } \min(u_{i,K}^p, u_{i,K\sigma}^p) \le 0, \\ \frac{u_{i,K}^p - u_{i,K\sigma}^p}{\log(u_{i,K}^p) - \log(u_{i,K\sigma}^p)} & \text{if } 0 \le u_{i,K}^p = u_{i,K\sigma}^p, \\ \frac{u_{i,K}^p - u_{i,K\sigma}^p}{\log(u_{i,K}^p) - \log(u_{i,K\sigma}^p)} & \text{otherwise.} \end{cases} $$

We also denote by $u_{\sigma}^p := (u_{i,\sigma}^p)_{1 \le i \le n}$. This choice for the edge concentration is crucial for the preservation at the discrete level of a discrete entropy-entropy dissipation inequality similar to (14) on the continuous level. Setting $u_{i,\sigma}^p = 0$ when $u_{i,K}^p$ or $u_{i,K\sigma}^p$ take negative values allows to extend in a continuous way the logarithmic mean to the whole $\mathbb{R}^2$. Moreover, this seemingly non-consistent choice allows to establish the positivity of the volume fractions by a simple contradiction argument, cf. Lemma 3.1 below. Therefore, this situation never occurs in practice, so that $u_{i,\sigma}^p$ is a genuine logarithmic mean.

The conservation laws are discretized in a conservative way with a time discretization relying on the backward Euler scheme:

$$ (34b) \qquad m_K \frac{u_{i,K}^p - u_{i,K}^{p-1}}{\Delta t_p} + \sum_{\sigma \in \mathcal{E}_K} m_{\sigma} J_{i,K\sigma}^p = 0, \quad \forall K \in \mathcal{T}, \forall 1 \le i \le n. $$

The relation between the fluxes and the variations of the volume fractions across the edges relies on formula (21) rather than on (5). This trick takes its inspiration in [18], and appears to be crucial in what follows for the derivation of the discrete
counterpart of the entropy-entropy dissipation estimate (14). More precisely, the discrete fluxes $J_{K\sigma}^p := (J_{i,K\sigma}^p)_{1\le i\le n}$ are solution to the following set of equations: for all $K \in \mathcal{T}$ and $\sigma \in \mathcal{E}_{K,\text{int}}$,

$$ \frac{1}{d_{\sigma}} D_{K\sigma} \mathbf{u}_i^p + c^* J_{i,K\sigma}^p + \sum_{1 \le j \le n} \bar{A}_{ij}(u_{\sigma}^p) J_{j,K\sigma}^p = 0, \quad \forall 1 \le i \le n, $$

which rewrites in a more compact form as

$$ (34c) \qquad \frac{1}{d_{\sigma}} D_{K\sigma} \mathbf{u}^p + c^* J_{K\sigma}^p + \bar{A}(u_{\sigma}^p) J_{K\sigma}^p = 0. $$

The discrete counterpart to the no-flux boundary condition (7) is naturally

$$ (34d) \qquad J_{K\sigma}^p = 0, \quad \forall \sigma \in \mathcal{E}_{K,\text{ext}}, K \in \mathcal{T}, 1 \le p \le P_T. $$

Notice that Formula (34c) yields conservative fluxes, i.e.,

$$ (35) \qquad J_{K\sigma}^p + J_{L\sigma}^p = 0, \quad \forall \sigma = K|L \in \mathcal{E}_{\text{int}}, 1 \le p \le P_T. $$

**Remark 2.2.** We stress on the fact here that we do not impose the constraint $J_{K\sigma}^p \in \mathcal{V}_0$ for all $K \in \mathcal{T}$, $\sigma \in \mathcal{E}_K$, and $1 \le p \le P_T$. Indeed, (34c) can be rewritten equivalently as

$$ \frac{1}{d_{\sigma}} D_{K\sigma} \mathbf{u}^p + (c^* I + \bar{A}(u_{\sigma}^p)) J_{K\sigma}^p = 0, $$

and the matrix $c^* I + \bar{A}(u_{\sigma}^p)$ differs in general from $A(u_{\sigma}^p)$ since $u_{\sigma}^p$ does not belong to $\mathcal{V}_1$ in general. As a consequence, Ker $(c^* I + \bar{A}(u_{\sigma}^p))$ may not be of dimension 1. Actually, we will see in Lemma 3.1 and Lemma 3.2 that for any $\mathbf{u}^{p-1}$ in

$$ \mathcal{A}^{\mathcal{T}} = \left\{ \mathbf{v} \in (\mathbb{R}_{+}^{\mathcal{T}})^n \mid (v_{i,K})_{1 \le i \le n} \in \mathcal{A} \text{ for all } K \in \mathcal{T} \right\}, $$

then any solution $\mathbf{u}^p$ to the scheme presented above also belongs to $\mathcal{A}^{\mathcal{T}}$ and that there exists a unique set of fluxes $(J_{K\sigma}^p)_{K \in \mathcal{T}, \sigma \in \mathcal{E}_K}$ satisfying (34c)-(34d), and that $J_{K\sigma}^p$ necessarily belongs to $\mathcal{V}_0$.

## 2.3. Main results and organisation.
We gather the main results of our paper in this section. Our first theorem concerns the existence of a discrete solution for a given mesh, and the preservation of the structural properties listed in Section 1.2.

In order to obtain a discrete counterpart of the entropy-entropy dissipation inequality (14), we need to introduce the discrete entropy functional $E_{\mathcal{T}} : (\mathbb{R}_+^{\mathcal{T}})^n \to \mathbb{R}$, which is defined by

$$ (36) \qquad E_{\mathcal{T}}(\mathbf{v}) = \sum_{i=1}^{n} \sum_{K \in \mathcal{T}} m_K v_{i,K} \log(v_{i,K}), \quad \forall \mathbf{v} = (v_i)_{1 \le i \le n} \in (\mathbb{R}_+^{\mathcal{T}})^n. $$

Note that the functional $E_{\mathcal{T}}$ is uniformly bounded on the set $\mathcal{A}^{\mathcal{T}}$. More precisely, there holds

$$ (37) \qquad -m_{\Omega} \log(n) \le E_{\mathcal{T}}(\mathbf{v}) \le 0, \quad \forall \mathbf{v} \in \mathcal{A}^{\mathcal{T}}. $$

Denote by $\mathbf{1}_{\mathcal{T}} = (1, \dots, 1) \in \mathbb{R}^{\mathcal{T}}$, then the following theorem holds:

**Theorem 2.3.** Let $(\mathcal{T}, \mathcal{E}, (x_K)_{K \in \mathcal{T}})$ be an admissible mesh and let $\mathbf{u}^0$ be defined by (33) from an initial condition $u^0 \in L^{\infty}(\Omega; \mathcal{A})$ satisfying the nondegeneracy assumption (9). Then, for all $1 \le p \le P_T$, the nonlinear system of equations (34) has (at least) a (strictly) positive solution $\mathbf{u}^p \in \mathcal{A}^{\mathcal{T}}$. This solution $\mathbf{u}^p$ satisfies
$\langle \mathbf{u}^p, \mathbf{1}_T \rangle_T = M$ and the corresponding fluxes $\mathbf{J}^p = (J_{K\sigma}^p)_{\sigma \in \mathcal{E}}$ are uniquely determined by (34c)-(34d) and belong to $(\mathcal{V}_0)^\mathcal{E}$, i.e. $\sum_{i=1}^n J_{i,K\sigma}^p = 0$ for all $\sigma \in \mathcal{E}$. Moreover, the following entropy-entropy dissipation estimate holds:

$$ (38) \quad E_T(\mathbf{u}^p) + \Delta t_p \sum_{\sigma=K|L \in \mathcal{E}_{\text{int}}} \left( \frac{c^*}{2} m_\sigma d_\sigma |J_{K\sigma}^p|^2 + \frac{\alpha}{2} \tau_\sigma |D_{K\sigma} \sqrt{\mathbf{u}^p}|^2 \right) \le E_T(\mathbf{u}^{p-1}). $$

The proof of Theorem 2.3 will be the purpose of Section 3.

From an iterated discrete solution $(\mathbf{u}, \mathbf{J}) = (\mathbf{u}^p, \mathbf{J}^p)_{1\le p \le P_T}$ to the scheme (34), we define for all $1 \le i \le n$, the piecewise constant approximate volume fractions $u_{i,\mathcal{T},\Delta t} : Q_T \to (0,1)$ defined almost everywhere by

$$ (39) \quad u_{i,\mathcal{T},\Delta t}(t,x) = u_{i,K}^n \quad \text{if } (t,x) \in (t_{p-1}, t_p] \times K. $$

Since $\mathbf{u}^p \in \mathcal{A}^T$, then $u_{\mathcal{T},\Delta t} = (u_{i,\mathcal{T},\Delta t})_{1\le i \le n}$ belongs to $L^\infty(Q_T; \mathcal{A})$. We also define approximate fluxes $J_{\mathcal{E},\Delta t} = (J_{i,\mathcal{E},\Delta t})_{1\le i \le n} : Q_T \to (\mathcal{V}_0)^d$ from the discrete fluxes $\mathbf{J}^p$ by setting

$$ (40) \quad J_{\mathcal{E},\Delta t}(t,x) = d J_{K\sigma}^p \nu_{K\sigma} \quad \text{if } (t,x) \in (t_{p-1}, t_p] \times \Delta_\sigma. $$

We are now in position to present our second main result, which concerns the convergence of the scheme as the discretisation parameters tend to 0. In what follows, let $(\mathcal{T}_m, \mathcal{E}_m, (x_K)_{K \in \mathcal{T}_m})_{m \ge 1}$ and $(\Delta \mathbf{t}_m)_{m \ge 1}$ be sequences of admissible discretisations of $\Omega$ and $(0, T)$ respectively. We assume that

$$ (41) \quad h_{\mathcal{T}_m} \underset{m \to \infty}{\longrightarrow} 0, \quad h_{T_m} \underset{m \to \infty}{\longrightarrow} 0, \quad \text{while } \inf_{m \ge 1} \zeta_{\mathcal{T}_m} = \zeta^* > 0. $$

Then, the following theorem holds:

**Theorem 2.4.** Let $(\mathcal{T}_m, \mathcal{E}_m, (x_K)_{K \in \mathcal{T}_m})_{m \ge 1}$ and $(\Delta \mathbf{t}_m)_{m \ge 1}$ be sequences of admissible discretisations of $\Omega$ and $(0, T)$ respectively fulfilling condition (41). Let $(\mathbf{u}_m, \mathbf{J}_m)_m = ((\mathbf{u}^p, \mathbf{J}^p)_{1\le p \le P_T}, m)_{m \ge 1}$ be a corresponding sequence of discrete solutions to (34), from which a sequence of approximate solutions $(u_{\mathcal{T}_m}, \Delta \mathbf{t}_m, J_{\mathcal{E}_m, \Delta \mathbf{t}_m})_{m \ge 1}$ is reconstructed thanks to (39)-(40). Then there exists a weak solution $(u, J)$ to (2)-(5)-(6) in the sense of Definition 1.1 such that, up to a upsequence,

$$ u_{\mathcal{T}_m, \Delta \mathbf{t}_m} \underset{m \to +\infty}{\longrightarrow} u \quad \text{a.e. in } Q_T, $$

and

$$ J_{\mathcal{E}_m, \Delta \mathbf{t}_m} \underset{m \to +\infty}{\longrightarrow} J \quad \text{weakly in } L^2((0, T) \times \Omega)^{d \times n}. $$

The proof of Theorem 2.4 is the purpose of Section 4. It is based on compactness arguments that are deduced from the a priori estimates established in Theorem 2.3.

Our third and last main result is about the long time behavior of the scheme. Since the finite time horizon $T$ is arbitrary in the above statement, we can consider the situation when $T = +\infty$ and study the long time behavior of the scheme. At the continuous level, it is shown in [41] that

$$ \|u(t, \cdot) - u^\infty\|_{L^1(\Omega)^n} \le C e^{-\lambda t}, \quad \forall t > 0, $$

where $u^\infty = (u_i^\infty)_{1\le i \le n}$ is the constant in space state with same mass as $u^0$, that is $u_i^\infty = M_i/m_\Omega$, $\bar{M}_i$ being defined by (9), where $C$ only depends on $\Omega$ and of the number of species $n$, and where $\lambda$ only depends on $\Omega$ and on $\alpha$ defined by (31). This result has a the following discrete counterpart.
**Theorem 2.5.** There exists C only depending on $\Omega$ and of the number of species $n$, and $\lambda$ only depending on $\Omega$ and on $\alpha$ such that

$$ (42) \qquad \|u_{\mathcal{T},\Delta t}(t, \cdot) - u^{\infty}\|_{L^1(\Omega)^n} \le Ce^{-\lambda t}, \quad \forall t > 0, $$

The proof of Theorem 2.5 is not difficult once all the material for the proofs of Theorems 2.3 and 2.4 has been introduced. We provide in Section 5.2 some numerical evidence for Theorem 2.5 as well as a sketch of its proof.

### 3. NUMERICAL ANALYSIS ON A FIXED GRID

This section is devoted to the proof of Theorem 2.3.

**3.1. A priori estimates.** The first lemma shows the non-negativity and the mass conservation of the solution to (34), together with the uniqueness of associated fluxes.

**Lemma 3.1.** Given $u^{p-1} \in \mathcal{A}^T$ satisfying

$$ (43) \qquad \langle \mathbf{1}_{\mathcal{T}}, \mathbf{u}^{p-1} \rangle_{\mathcal{T}} = M \in (\mathbb{R}_+^*)^n, $$

then any solution $u^p$ to (34) satisfies $\langle \mathbf{1}_{\mathcal{T}}, \mathbf{u}^p \rangle_{\mathcal{T}} = M$ and is positive in the sense that $u_{i,K}^p > 0$ for all $K \in \mathcal{T}$ and all $1 \le i \le n$. Besides, for any solution $u^p$ to (34), there exists a unique set of fluxes $J^p$ satisfying (34c)-(34d).

*Proof.* Let $u^p$ be a solution to (34) and let $1 \le i \le n$. Let us first prove that the total volume of each species is conserved, so that $\langle \mathbf{1}_{\mathcal{T}}, \mathbf{u}^p \rangle_{\mathcal{T}} = M$. Summing equation (34b) over $K \in \mathcal{T}$ gives

$$ \langle \mathbf{1}_{\mathcal{T}}, \mathbf{u}^p \rangle_{\mathcal{T}} - \langle \mathbf{1}_{\mathcal{T}}, \mathbf{u}^{p-1} \rangle_{\mathcal{T}} = -\Delta t_p \sum_{\sigma=K|L \in \mathcal{E}_{\text{int}}} m_{\sigma} (J_{K\sigma}^p + J_{L\sigma}^p) - \Delta t_p \sum_{\sigma \in \mathcal{E}_{\text{ext}}} m_{\sigma} J_{K\sigma}^p. $$

Then it follows directly from the local conservativity of the scheme (35) and from the discrete no-flux boundary condition (34d) that

$$ \langle \mathbf{1}_{\mathcal{T}}, \mathbf{u}^p \rangle_{\mathcal{T}} = \langle \mathbf{1}_{\mathcal{T}}, \mathbf{u}^{p-1} \rangle_{\mathcal{T}} = M. $$

Let us now prove that $u^p$ is positive. Let $1 \le i \le n$. We consider a cell $K \in \mathcal{T}$ where $u_i^p$ reaches its minimum, i.e. such that $u_{i,K}^p \le u_{i,L}^p$ for all $L \in \mathcal{T}$. Assume for contradiction that $u_{i,K}^p \le 0$. Let us recall again equation (34b), which implies that

$$ (44) \qquad m_K \frac{u_{i,K}^p - u_{i,K}^{p-1}}{\Delta t_p} = - \sum_{\sigma \in \mathcal{E}_K} m_{\sigma} J_{i,K\sigma}^p. $$

On the one hand, the term on the left-hand side is non-positive since $u_{i,K}^{p-1} \ge 0 \ge u_{i,K}^p$. On the other hand, the specific choice (34a) for the edge volume fractions implies that $u_{i,\sigma}^p = 0$ for all $\sigma \in \mathcal{E}_K$. Therefore, $\bar{A}_{ij}(u_{\sigma}^p) = 0$ for all $1 \le j \ne i \le n$. As a consequence, relation (34c) reduces to

$$ \frac{1}{d_{\sigma}} D_{K\sigma} \mathbf{u}_{i}^{p} + \left( c^{*} + \sum_{1 \le j \ne i \le n} \bar{c}_{ij} u_{j,\sigma}^{p} \right) J_{i,K\sigma}^{p} = 0. $$

Since $D_{K\sigma}u_i^p \ge 0$ and $u_{j,\sigma}^p \ge 0$, $\bar{c}_{ij} \ge 0$ for all $1 \le j \le n$, we obtain that $J_{i,K\sigma}^p \le 0$ for all $\sigma \in \mathcal{E}_K$. Using (44), this yields that $J_{i,K\sigma}^p = \tau_{\sigma} D_{K\sigma}u_i^p = 0$ for all $\sigma \in \mathcal{E}_K$.
As a consequence, $u_{i,K}^p = u_{i,L}^p$ for all $L \in \mathcal{T}$ such that $\sigma = K|L \in \mathcal{E}_K$. Iterating this argument and since $\Omega$ is connected, we thus obtain that $u_{i,L}^p = u_{i,K}^p \le 0$ for all $L \in \mathcal{T}$. This implies that $\langle \mathbf{u}_i^p, \mathbf{1}_\mathcal{T} \rangle_\mathcal{T} \le 0$ which contradicts the property $\langle \mathbf{1}_\mathcal{T}, \mathbf{u}_i^p \rangle_\mathcal{T} = M_i > 0$ we just established. Thus, $\mathbf{u}^p$ is positive.

As a consequence, for all $\sigma \in \mathcal{E}_{\text{int}}$ and all $1 \le i \le n$, $u_{i,\sigma}^p > 0$. The fact that there exists a unique $\mathbf{J}^p$ associated to $\mathbf{u}^p$ via (34c)-(34d) is then a consequence of Lemma 1.2. Indeed, for all $K \in \mathcal{T}$ and all $\sigma \in \mathcal{E}_{K,\text{int}}$, noticing that $D_{K\sigma}\mathbf{u}^p = M(u_\sigma^p)D_{K\sigma}\log(\mathbf{u}^p)$ due to the specific choice (34a) for $u_\sigma^p$, we can rewrite equivalently (34c) as

$$ \frac{1}{d_{\sigma}} M(u_{\sigma}^{p}) D_{K\sigma} \log(\mathbf{u}^{p}) + (c^{*}I + \bar{A}(u_{\sigma}^{p})) J_{K\sigma}^{p} = 0. $$

The positivity of $u_\sigma^p$ implies the inversibility of matrix $M(u_\sigma^p)$. As a consequence, it holds that

$$ (45) \qquad \frac{1}{d_{\sigma}} D_{K\sigma} \log(\mathbf{u}^p) + (c^* M(u_{\sigma}^p)^{-1} + M(u_{\sigma}^p)^{-1} \bar{A}(u_{\sigma}^p)) J_{K\sigma}^p = 0. $$

Moreover, thanks to Lemma 1.2, the matrix $B(u_\sigma^p) = c^*M(u_\sigma^p)^{-1} + M(u_\sigma^p)^{-1}\bar{A}(u_\sigma^p)$ is symmetric positive definite, and the only solution $J_{K\sigma}^p$ to (45) is given by

$$ (46) \qquad J_{K\sigma}^p = -\frac{1}{d_{\sigma}} B(u_{\sigma}^p)^{-1} D_{K\sigma} \log(\mathbf{u}^p). $$

Hence the desired result.

□

The next lemma shows that the total discrete flux vanishes across all edges and that the volume filling constraint is automatically satisfied without being enforced.

**Lemma 3.2.** Given $\mathbf{u}^{p-1} \in \mathcal{A}^T$ satisfying (43), any solution $(\mathbf{u}^p, \mathbf{J}^p)$ to (34) belongs to $\mathcal{A}^T \times (V_0)^\mathcal{E}$.

*Proof.* Since $\mathbf{u}^{p-1} \in \mathcal{A}^T$ satisfies (43), $\mathbf{u}^{p-1}$ is nonnegative, and using Lemma 3.1, any corresponding solution $\mathbf{u}^p$ to (34) is then positive.

Let us define $\mathbf{w}^{p-1} = (w_K^{p-1})_{K \in \mathcal{T}} := \langle \mathbf{1}, \mathbf{u}^{p-1} \rangle$ and analogously $\mathbf{w}^p = \langle \mathbf{1}, \mathbf{u}^p \rangle$, and let us denote by $G_{\sigma K}^p := \langle \mathbf{1}, J_{K\sigma}^p \rangle$ for all $K \in \mathcal{T}$ and $\sigma \in \mathcal{E}_K$. Summing equations (34b) for $i = 1, \dots, n$, we obtain that

$$ m_K \frac{w_K^p - w_K^{p-1}}{\Delta t_p} = - \sum_{\sigma \in \mathcal{E}_K} m_\sigma G_{K\sigma}^p. $$

In addition, summing (34c) over $i$ provides that for all $\sigma = K|L \in \mathcal{E}_{\text{int}}$,

$$ 0 = \frac{1}{d_{\sigma}} D_{K\sigma} \mathbf{w}^p + c^* G_{K\sigma}^p + \langle \mathbf{1}, \bar{A}(u_{\sigma}^p) J_{K\sigma}^p \rangle \stackrel{(19)}{=} \frac{1}{d_{\sigma}} D_{K\sigma} \mathbf{w}^p + c^* G_{K\sigma}^p. $$

Thus, $\mathbf{w}$ is solution to the classical backward Euler TPFA scheme for the heat equation with diffusion coefficient $\frac{1}{c^*}$. This scheme is well-posed and $\mathbf{w}^p = \mathbf{w}^{p-1} = \mathbf{1}_\mathcal{T}$ is its unique solution, which implies that $\mathbf{u} \in \mathcal{A}^T$. Moreover, the fluxes $G_{K\sigma}^p$ are all equal to zero, so that $J_{K\sigma} \in (V_0)^\mathcal{E}$. □

**Remark 3.3.** A direct consequence of Lemma 3.1 is that the edge volume fractions originally defined by (34a) are genuine logarithmic means of the surrounding cell volume fractions. Since the logarithmic mean is smaller than the arithmetic one, one readily checks that

$$ (47) \qquad \langle \mathbf{1}, u_{\sigma}^{p} \rangle \le 1, \quad \forall \sigma \in \mathcal{E}_{\text{int}}, \forall p \ge 1. $$
The last statement of this section is devoted to the entropy entropy-dissipation estimate (38).

**Lemma 3.4.** Given $\mathbf{u}^{p-1} \in \mathcal{A}^{\mathcal{T}}$, any solution $(\mathbf{u}^p, \mathbf{J}^p) \in \mathcal{A}^{\mathcal{T}} \times (\mathcal{V}_0)^{\mathcal{E}}$ to (34) satisfies

$$ (48) \quad E_{\mathcal{T}}(\mathbf{u}^p) + \Delta t_p \sum_{\sigma=K|L \in \mathcal{E}_{\text{int}}} \left( \frac{c^*}{2} m_{\sigma} d_{\sigma} |J_{K\sigma}^p|^2 + \frac{\alpha}{2} \tau_{\sigma} |D_{K\sigma} \sqrt{\mathbf{u}^p}|^2 \right) \le E_{\mathcal{T}}(\mathbf{u}^{p-1}). $$

*Proof.* Multiplying equation (34b) by $\Delta t_p \log(u_{i,K}^p)$ (which makes sense since $\mathbf{u}^p$ is positive owing to Lemma 3.1), and summing over all the cells and species leads to

$$ (49) \qquad T_1 + T_2 = 0, $$

where we have set

$$ T_1 = \sum_{K \in \mathcal{T}} \sum_{i=1}^{n} \left[ u_{i,K}^{p} \log(u_{i,K}^{p}) - u_{i,K}^{p-1} \log(u_{i,K}^{p}) \right] m_K, $$

$$ T_2 = \Delta t_p \sum_{i=1}^{n} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} m_{\sigma} J_{i,K\sigma}^{p} \log(u_{i,K}^{p}). $$

On the one hand, using the convexity of the function $\mathbb{R}_+ \ni x \mapsto x \log x$, it holds that

$$ u_{i,K}^{p} - u_{i,K}^{p-1} + u_{i,K}^{p} \log(u_{i,K}^{p}) - u_{i,K}^{p-1} \log(u_{i,K}^{p}) \geq u_{i,K}^{p} \log(u_{i,K}^{p}) - u_{i,K}^{p-1} \log(u_{i,K}^{p-1}), $$

which implies, together with Lemma 3.2, that

$$ (50) \qquad T_1 \geq E_{\mathcal{T}}(\mathbf{u}^p) - E_{\mathcal{T}}(\mathbf{u}^{p-1}). $$

On the other hand, the conservativity of the fluxes (35) and the discrete no-flux boundary condition (34d) allow to reorganise the term $T_2$ as

$$ T_2 = -\Delta t_p \sum_{\sigma=K|L \in \mathcal{E}_{\text{int}}} m_{\sigma} \langle J_{K\sigma}^p, D_{K\sigma} \log(\mathbf{u}^p) \rangle. $$

Bearing in mind the expression (46) of the fluxes,

$$ - \langle J_{K\sigma}^p, D_{K\sigma} \log(\mathbf{u}^p) \rangle = \frac{d_{\sigma}}{2} \langle J_{K\sigma}^p, B(u_{\sigma}^p) J_{K\sigma}^p \rangle \\ + \frac{1}{2d_{\sigma}} \langle D_{K\sigma} \log(\mathbf{u}^p), B(u_{\sigma}^p)^{-1} D_{K\sigma} \log(\mathbf{u}^p) \rangle. $$

Then estimates (29) provide that

$$ \langle J_{K\sigma}^p, B(u_{\sigma}^p) J_{K\sigma}^p \rangle \geq c^* |J_{K\sigma}|^2 $$

and, recalling the definition (31) of $\alpha$,

$$ \langle D_{K\sigma} \log(\mathbf{u}^p), B(u_{\sigma}^p)^{-1} D_{K\sigma} \log(\mathbf{u}^p) \rangle \geq \frac{\alpha}{4} \langle D_{K\sigma} \log(\mathbf{u}^p), M(u_{\sigma}^p) D_{K\sigma} \log(\mathbf{u}^p) \rangle. $$

Thanks to the particular choice (34a) for $u_{\sigma}^p$, the right-hand side rewrites

$$ \langle D_{K\sigma} \log(\mathbf{u}^p), M(u_{\sigma}^p) D_{K\sigma} \log(\mathbf{u}^p) \rangle = \langle D_{K\sigma} \log(\mathbf{u}^p), D_{K\sigma} \mathbf{u}^p \rangle \geq 4 |D_{K\sigma} \sqrt{\mathbf{u}^p}|^2, $$

the last inequality being a consequence of the elementary inequality

$$ (a-b)(\log(a) - \log(b)) \geq 4(\sqrt{a} - \sqrt{b})^2 $$
holding for any positive a, b. Summing up, we have

$$ (51) \qquad T_2 \ge \Delta t_p \sum_{\sigma=K|L \in \mathcal{E}_{\text{int}}} \left( \frac{c^*}{2} m_\sigma d_\sigma |J_{K\sigma}^p|^2 + \frac{\alpha}{2} \tau_\sigma |D_{K\sigma} \sqrt{\mathbf{u}^p}|^2 \right). $$

To conclude the proof, it only remains to incorporate (50) and (51) in (49). □

## 3.2. Existence of discrete solutions.
The purpose of this section is to prove the existence of a solution to (34).

**Proposition 3.5.** Given $\mathbf{u}^{p-1} \in \mathcal{A}^T$ satisfying (43), then there exists at least one solution $(\mathbf{u}^p, \mathbf{J}^p) \in \mathcal{A}^T \times (\mathcal{V}_0)^\mathcal{E}$ to the scheme (34).

*Proof.* The proof relies on a topological degree argument [44, 24]. The idea is to transform continuously our complex nonlinear system into a linear system while guaranteeing that enough a priori estimates controlling the solution remain valid all along the homotopy. We sketch the main ideas of the proof, making the homotopy explicit.

For $\lambda \in [0, 1]$, we look for $(\mathbf{u}^{(\lambda)}, \mathbf{J}^{(\lambda)}) \in \mathbb{R}^{n \times T} \times \mathbb{R}^{n \times \mathcal{E}}$ solution to the algebraic system (34) where the matrix $\bar{A}(u_\sigma^{(\lambda)})$ is replaced by $\lambda \bar{A}(u_\sigma^{(\lambda)})$. Our system (34) corresponds to the case $\lambda = 1$, whereas the case $\lambda = 0$ corresponds to the usual TPFA finite volume scheme for $n$ decoupled heat equations all with the same diffusion coefficient $\frac{1}{c^*}$. Mimicking the calculations presented in Section 3.1, one shows that whatever $\lambda \in [0, 1]$, any corresponding solution $(\mathbf{u}^{(\lambda)}, \mathbf{J}^{(\lambda)})$ lies in $\mathcal{A}^T \times (\mathcal{V}_0)^\mathcal{E}$, and $\mathbf{u}^{(\lambda)}$ is positive. Moreover, the entropy - entropy dissipation estimate and the uniform bound (37) on the entropy ensure that

$$ \|\mathbf{J}^{(\lambda)}\|_{\mathcal{E}}^2 \le \frac{2m_{\Omega} \log n}{c^* \Delta t_p} =: K. $$

where $\|\mathbf{J}^{(\lambda)}\|_{\mathcal{E}}^2 = \sum_{\sigma=K|L \in \mathcal{E}_{\text{int}}} m_\sigma d_\sigma |J_{K\sigma}^{(\lambda)}|^2$. Fixing $\eta > 0$, we define the relatively compact open sets

$$ \mathcal{A}_{\eta}^{T} = \left\{ \mathbf{u} \in (\mathbb{R}^{T})^{n} \mid \inf_{\mathbf{v} \in \mathcal{A}^{T}} \|\mathbf{u} - \mathbf{v}\| < \eta \right\} $$

and

$$ (\mathcal{V}_0)_\eta^\mathcal{E} = \left\{ \mathbf{J} \in (\mathbb{R}^\mathcal{E})^n \mid \|\mathbf{J}\|_{\mathcal{E}}^2 < K^{1/2} + \eta \quad \text{and} \quad \inf_{\mathbf{F} \in (\mathcal{V}_0)^\mathcal{E}} \|\mathbf{J} - \mathbf{F}\| < \eta \right\}. $$

The a priori estimates ensure that no solution $(\mathbf{u}^{(\lambda)}, \mathbf{J}^{(\lambda)})$ of the modified scheme can cross the boundary of the open set $\mathcal{A}_{\eta}^T \times (\mathcal{V}_0)_\eta^\mathcal{E}$. The topological degree associated to the modified scheme and $\mathcal{A}_{\eta}^T \times (\mathcal{V}_0)_\eta^\mathcal{E}$ is constant with respect to $\lambda$, and takes the value +1 for $\lambda = 0$ since the system is linear and invertible with positive determinant. So it is also equal to 1 for $\lambda = 1$, ensuring the existence of a solution to the nonlinear problem (34). □

The proof of Theorem 2.3 is now complete.
## 4. PROOF OF THEOREM 2.4

We consider here a sequence $(\mathcal{T}_m, \mathcal{E}_m, (x_K)_{K \in \mathcal{T}_m})_{m \ge 1}$ of admissible space discretizations with $h_{\mathcal{T}_m}$ going to 0 as $m$ tends to $+\infty$, while the regularity $\zeta_{\mathcal{T}_m}$ remains uniformly bounded from below by a positive constant $\zeta^*$. We also consider a sequence $(\Delta \mathbf{t}_m)_{m \ge 1} = ((\Delta t_{p,m})_{1 \le p \le P_{T,m}})_{m \ge 1}$ of admissible time discretizations such that $h_{T,m}$ goes to 0 as $m$ goes to infinity.

From the discrete solutions $(\mathbf{u}_m, \mathbf{J}_m)$, $m \ge 1$, the existence of which being guaranteed by Theorem 2.3, we reconstruct the piecewise constant functions $u_{\mathcal{T}_m, \Delta \mathbf{t}_m} \in L^\infty(Q_T; A)$ and $J_{\mathcal{E}_m, \Delta \mathbf{t}_m} \in L^2(Q_T; V_0)^d$ thanks to formulas (39) and (40). In the convergence analysis, we also need the weakly consistent piecewise constant gradient reconstruction operators $\nabla_{\mathcal{E}_m}$ and $\nabla_{\mathcal{E}_m, \Delta \mathbf{t}_m}$ defined for $m \ge 1$ and $\mathbf{v} \in \mathbb{R}^{\mathcal{T}_m}$

$$ (52) \qquad \nabla_{\mathcal{E}_m} \mathbf{v}(x) = d D_{K_\sigma} \mathbf{v}_m \nu_{K_\sigma} \quad \text{if } x \in \Delta_\sigma, \sigma \in \mathcal{E}_m, $$

and, for $\mathbf{v} = (\mathbf{v}^p)_{0 \le p \le P_{T,m}} \in \mathbb{R}^{(1+P_{T,m}) \times \mathcal{T}_m}$,

$$ (53) \qquad \nabla_{\mathcal{E}_m, \Delta \mathbf{t}_m} \mathbf{v}_m(t, \cdot) = \nabla_{\mathcal{E}_m} \mathbf{v}^p \quad \text{if } t \in (t_{p-1}, t_p], \quad 1 \le p \le P_{T,m}. $$

### 4.1. Compactness on approximate reconstructions.
The next proposition is the main result of this section.

**Proposition 4.1.** There exists $u \in L^\infty(Q_T; A^T) \cap L^2(0, T; H^1(\Omega))^n$ with $\sqrt{u} \in L^2(0, T; H^1(\Omega))^n$, and $J \in L^2(Q_T; (V_0)^d)$ such that, up to a subsequence, the following convergence properties hold:

$$ (54) \qquad u_{\mathcal{T}_m, \Delta \mathbf{t}_m} \underset{m \to +\infty}{\longrightarrow} u \quad \text{a.e. in } Q_T, $$

$$ (55) \qquad \nabla_{\mathcal{E}_m, \Delta \mathbf{t}_m} \sqrt{\mathbf{u}_m} \underset{m \to +\infty}{\longrightarrow} \nabla \sqrt{u} \quad \text{weakly in } L^2(Q_T)^{n \times d}, $$

$$ (56) \qquad \nabla_{\mathcal{E}_m, \Delta \mathbf{t}_m} \mathbf{u}_m \underset{m \to +\infty}{\longrightarrow} \nabla u \quad \text{weakly in } L^2(Q_T)^{n \times d}, $$

$$ (57) \qquad J_{\mathcal{E}_m, \Delta \mathbf{t}_m} \underset{m \to +\infty}{\longrightarrow} J \quad \text{weakly in } L^2(Q_T)^{n \times d}. $$

*Proof.* Summing (38) over $p \in \{1, \dots, P_{T,m}\}$ and using the bound (37) on $E_T$ provides

$$ (58) \qquad \sum_{p=1}^{P_{T,m}} \Delta t_p \sum_{\sigma \in \mathcal{E}_{\text{int},m}} \left( \frac{\alpha}{2} \tau_{\sigma} |D_{\sigma} \sqrt{\mathbf{u}_m}|^2 + \frac{c^*}{2} m_{\sigma} d_{\sigma} |J_{K_{\sigma}}|^2 \right) \le m_{\Omega} \log n. $$

Recalling the elementary geometrical relation $dm_{\Delta_{\sigma}} = m_{\sigma} d_{\sigma}$ and the definitions (40) of $J_{\mathcal{E}_m, \Delta \mathbf{t}_m}$ and (52)-(53), one obtains that

$$ (59) \qquad \|J_{\mathcal{E}_m, \Delta \mathbf{t}_m}\|_{L^2(Q_T)^{n \times d}} + \|\nabla_{\mathcal{E}_m, \Delta \mathbf{t}_m} \sqrt{\mathbf{u}_m}\|_{L^2(Q_T)^{n \times d}} \le C $$

for some $C$ not depending on $m$. As a straightforward consequence, there exists $J, F \in L^2(Q_T)^{n \times d}$ such that (57) holds, as well as

$$ (60) \qquad \nabla_{\mathcal{E}_m, \Delta \mathbf{t}_m} \sqrt{\mathbf{u}_m} \underset{m \to +\infty}{\longrightarrow} F \quad \text{weakly in } L^2(Q_T)^{n \times d}. $$

The fact that $J \in L^2(Q_T; V_0)^d$ results from the stability of linear space $V_0$ for the weak convergence. Moreover, since $0 \le u_K^n \le 1$, then $D_{\sigma} \mathbf{u}_m^p \le 2D_{\sigma} \sqrt{\mathbf{u}_m^p}$ for all $\sigma \in \mathcal{E}_{\text{int},m}$ and all $1 \le p \le P_{T,m}$. Therefore, we deduce from (59) that

$$ \|\nabla_{\mathcal{E}_m, \Delta \mathbf{t}_m} \mathbf{u}_m\|_{L^2(Q_T)^{n \times d}} \le C, $$
whence the existence of some $G \in L^2(Q_T)$ such that

$$ (61) \qquad \nabla_{\mathcal{E}_m, \Delta t_m} u_m \underset{m \to +\infty}{\longrightarrow} G \quad \text{weakly in } L^2(Q_T)^{n \times d}. $$

On the other hand, $u_{\mathcal{T}_m, \Delta t_m}$ belongs to the bounded subset $L^\infty(Q_T; \mathcal{A})$ of $L^\infty(Q_T)^n$ for all $m \ge 1$. Therefore, up to a subsequence, $u_{\mathcal{T}_m, \Delta t_m}$ converges in the $L^\infty(Q_T)^n$-weak star sense towards some $u$, which takes its values in $\mathcal{A}$ since both the positivity and the sum to 1 property are stable when passing to the limit in this topology.

To conclude this proof, it remains to check that the convergence of $u_{\mathcal{T}_m, \Delta t_m}$ towards $u$ holds point-wise, and to identify $F$ and $G$ as $\nabla\sqrt{u}$ and $\nabla u$ respectively. These properties are provided all at once by the nonlinear discrete Aubin-Simon lemma [4, Theorem 3.9]. As already established in [4], this theorem applies naturally in the TPFA finite volume context. The only point to be checked is a discrete $L^2(0, T; H^{-1}(\Omega))$ estimate on the time increments of $u_{\mathcal{T}_m, \Delta t_m}$. More precisely, for $\phi \in C_c^\infty((0, T) \times \Omega; \mathbb{R}^n)$, one defines $\phi = (\phi_{i,K}^p) \in \mathbb{R}^{n \times P_{T,m} \times \mathcal{T}_m}$ by

$$ \phi_{i,K}^p = \frac{1}{\Delta t_p m_K} \int_{t_{p-1}}^{t_p} \int_K \phi_i(t, x) dx dt. $$

It follows from (34b)-(34d)-(35) that

$$ \sum_{p=1}^{P_{T,m}} \sum_{K \in \mathcal{T}_m} m_K \langle (u_K^p - u_K^{p-1}), \phi_K^p \rangle = \sum_{p=1}^{P_{T,m}} \Delta t_p \sum_{\sigma \in \mathcal{E}_{\text{int},m}} m_\sigma \langle J_{K\sigma}^p, D_{K\sigma} \phi^p \rangle. $$

Applying Cauchy-Schwarz inequality leads to

$$ \sum_{p=1}^{P_{T,m}} \sum_{K \in \mathcal{T}_m} m_K \langle (u_K^p - u_K^{p-1}), \phi_K^p \rangle \\ \le \left( \sum_{p=1}^{P_{T,m}} \Delta t_p \sum_{\sigma \in \mathcal{E}_{\text{int},m}} m_\sigma d_\sigma |J_{K\sigma}^p|^2 \right)^{1/2} \left( \sum_{p=1}^{P_{T,m}} \Delta t_p \sum_{\sigma \in \mathcal{E}_{\text{int},m}} \tau_\sigma |D_\sigma \phi^p|^2 \right)^{1/2}. $$

The discrete $L^2(Q_T)^d$ estimate on the fluxes (58) shows that the first term in the right-hand side is bounded, whereas the second term is the discrete $L^2(0, T; H^1(\Omega))$ semi-norm of $\phi$. A straightforward generalisation of [28, Lemma 9.4] shows that

$$ \sum_{p=1}^{P_{T,m}} \Delta t_p \sum_{\sigma \in \mathcal{E}_{\text{int},m}} \tau_\sigma |D_\sigma \phi^p|^2 \le C \|\nabla \phi\|_{L^2(Q_T)^d}^2 $$

for some $C$ only depending on the regularity factor $\zeta^*$. Therefore,

$$ \sum_{p=1}^{P_{T,m}} \sum_{K \in \mathcal{T}_m} m_K \langle (u_K^p - u_K^{p-1}), \phi_K^p \rangle \le C \|\nabla \phi\|_{L^2(Q_T)^d} \le C \|\nabla \phi\|_{L^\infty(Q_T)^d}, $$

which is exactly the condition required to apply [4, Theorem 3.9], which provides (54)-(55)-(56) all at once, concluding the proof of Proposition 4.1. $\square$

For all $m \ge 1$, we introduce the diamond cell based reconstruction $u_{\mathcal{E}_m, \Delta t_m}$ of the volume fractions defined by

$$ u_{\mathcal{E}_m, \Delta t_m}(t, x) = u_{\sigma}^p \quad \text{if } (t, x) \in (t_{p-1}, t_p] \times \Delta_{\sigma}, \sigma \in \mathcal{E}_m, 1 \le p \le P_{T,m}, $$
where the $u_{\sigma}^p$ are given by (34a). The following lemma shows that both reconstructions $u_{\mathcal{E}_m, \Delta t_m}$ and $u_{\mathcal{T}_m, \Delta t_m}$ share the same limit $u$. The proof is omitted there since it is similar to the one of [18, Lemma 4.4].

**Lemma 4.2.** Let $u$ be as in Proposition 4.1 then, up to a subsequence, $u_{\mathcal{E}_m, \Delta t_m}$ converges in $L^r(Q_T)$, $1 \le r < +\infty$ towards $u$ as $m$ tends to $+\infty$.

4.2. **Convergence towards a weak solution.** Our last statement to conclude the proof of Theorem 2.4 consists in identifying the limit values $(u, J)$ of the approximate solutions as weak solutions to the Stefan-Maxwell cross-diffusion system.

**Proposition 4.1.** Let $(u, J)$ be as in Proposition 4.1 then $(u, J)$ is a weak solution to (2)-(5)-(6) in the sense of Definition 1.1.

*Proof.* One has already established in Proposition 4.1 that the limit values $(u, J)$ lie in the right functional spaces. It only remains to check that (2), (7) and (21) hold in the distributional sense.

Equation (34c) implies that

$$ (62) \qquad \nabla_{\mathcal{E}_m, \Delta t_m} u_m + (c^* I + \bar{A}(u_{\mathcal{E}_m, \Delta t_m})) J_{\mathcal{E}_m, \Delta t_m} = 0, \quad \forall m \ge 1. $$

Since $v \mapsto \bar{A}(v)$ is continuous, it follows from Lemma 4.2 that $\bar{A}(u_{\mathcal{E}_m, \Delta t_m})$ tends to $A(u)$ in $L^2(Q_T)^{n \times n}$. Then thanks to the convergence properties (56)-(57), one can pass to the weak limit in (62) to recover that (21) holds in $L^1(Q_T)^{n \times d}$, thus also in $L^2(Q_T)^{n \times d}$.

Concerning equations (2) and (7), we establish them in the distributional sense (10). Let $\phi \in C_c^\infty([0, T) \times \bar{\Omega})$, then for $m \ge 1$, define $\phi_m = (\phi_K^p)_{K \in \mathcal{T}_m, 1 \le p \le P_{T,m}}$ by setting $\phi_K^p = \phi(t_p, x_K)$. Multiplying (2) by $\Delta t_p \phi_K^{p-1}$ for some $1 \le i \le N$ and summing over $K \in \mathcal{T}_m$ and $1 \le p \le P_{T,m}$ gives after reorganisation that

$$ (63) \qquad \iint_{Q_T} u_{i, \mathcal{T}_m, \Delta t_m} \partial_t \phi + \int_{\Omega} u_i^0 \phi(0, \cdot) + \iint_{Q_T} J_{i, \mathcal{E}_m, \Delta t_m} \cdot \nabla \phi \\ = R_{1,m}(\phi) + R_{2,m}(\phi) + R_{3,m}(\phi), $$

where we have set

$$ R_{1,m}(\phi) = \sum_{p=1}^{P_{T,m}} \sum_{K \in \mathcal{T}_m} m_K u_{i,K}^p \left( \phi_K^p - \phi_K^{p-1} - \frac{1}{m_K} \int_{t_{p-1}}^{t_p} \int_K \partial_t \phi \right), $$

$$ R_{2,m}(\phi) = \sum_{K \in \mathcal{T}_m} m_K u_{i,K}^0 \left( \phi_K^0 - \frac{1}{m_K} \int_K \phi(0, \cdot) \right), $$

$$ R_{3,m}(\phi) = \sum_{p=1}^{P_{T,m}} \Delta t_p \sum_{\sigma \in \mathcal{E}_m} m_{\sigma} d_{\sigma} J_{i,K\sigma}^p \left( \frac{1}{d_{\sigma}} D_{K\sigma} \phi_m^{p-1} - \frac{1}{m_{\Delta_{\sigma}} \Delta t_p} \int_{t_{p-1}}^{t_p} \nabla \phi \cdot \nu_{K\sigma} \right). $$

It follows from the regularity of $\phi$ that

$$ \left| \phi_K^p - \phi_K^{p-1} - \frac{1}{m_K} \int_{t_{p-1}}^{t_p} \int_K \partial_t \phi \right| \le C \Delta t_p (h_{\mathcal{T}_m} + h_{\mathcal{T}_m}), $$

so that, using that $0 \le u_{i,K}^p \le 1$, we obtain that

$$ (64) \qquad |R_{1,m}(\phi)| \le C(h_{\mathcal{T}_m} + h_{\mathcal{T}_m}) \underset{m \to +\infty}{\longrightarrow} 0. $$
Similarly, one shows that

$$ (65) \qquad |R_{2,m}(\phi)| \le C h_{\tau_m} \underset{m \to +\infty}{\longrightarrow} 0. $$

Finally, the orthogonality condition on the mesh, namely point (iii) of Definition (2.1), ensures that

$$ \left| \frac{1}{d_{\sigma}} D_{K\sigma} \phi_m^{p-1} - \frac{1}{m_{\Delta_{\sigma}} \Delta t_p} \int_{t_{p-1}}^{t_p} \nabla \phi \cdot \nu_{K\sigma} \right| \le C(h_{\tau_m} + h_{T_m}). $$

Therefore,

$$ |R_{3,m}(\phi)| \le C(h_{\tau_m} + h_{T_m}) \|J_{i,\varepsilon_m, \Delta t_m}\|_{L^1(Q_T)^d} \underset{m \to +\infty}{\longrightarrow} 0 $$

since $\|J_{i,\varepsilon_m, \Delta t_m}\|_{L^1(Q_T)^d}$ can be controlled thanks to the Cauchy-Schwarz inequality by $T^{1/2}m_{\Omega}^{1/2} \|J_{i,\varepsilon_m, \Delta t_m}\|_{L^2(Q_T)^d}$ which is bounded thanks to (59). Then in view of the convergence in $L^1(Q_T)$ of $u_{i,\tau_m, \Delta t_m}$ towards $u_i$ and of the weak convergence in $L^2(Q_T)^d$ of $J_{i,\varepsilon_m, \Delta t_m}$ towards $J_i$, one can pass to the limit in (63) to recover that

$$ \iint_{Q_T} u_i \partial_t \phi + \int_{\Omega} u_i^0 \phi(0, \cdot) + \iint_{Q_T} J_i \cdot \nabla \phi = 0. $$

The weak formulation (10) is then recovered by summing over $i$. $\square$

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

## 5. NUMERICAL RESULTS

The aim of this section is to collect some numerical results obtained with the numerical scheme presented in the preceding sections. The numerical scheme has been implemented using Julia and the different codes used to produce the numerical tests presented below can be found at [26] (10.5281/zenodo.3934286). The nonlinear system is solved thanks to a modified Newton algorithm with stopping criterion $\|u^{p,k+1} - u^{p,k}\|_{\ell^{\infty}} < 10^{-12}$ where the superscript $k$ refers to the iteration of the Newton method. The obtained solution, denoted by $u^{p-2/3}$ is then projected onto $\mathcal{A}$ by setting:

$$ (66) \qquad u^{p-1/3} = \max(u^{p-2/3}, 10^{-12}) \quad \text{then} \quad u_{i,K}^{p} = \frac{u_{i,K}^{p-1/3}}{\sum_{i=1}^{n} u_{i,K}^{p-1/3}}. $$

5.1. **Convergence under grid refinement.** We first present some numerical results obtained on a one-dimensional test case, in order to illustrate the rate of convergence of the method with respect to the spatial discretization parameter. Here, $\Omega = (0, 1)$, and we consider a system composed of three different species ($n=3$). Two different initial conditions $u_0$ are considered:

* a smooth initial profile defined for $x \in (0, 1)$ by

$$ (67) \qquad u_1^0(x) = u_2^0(x) = \frac{1}{4} + \frac{1}{4} \cos(\pi x); $$

* a non-smooth initial profile defined for $x \in (0, 1)$ by

$$ (68) \qquad u_1^0(x) = \mathbb{1}_{[3/8,5/8]}(x), \quad u_2^0(x) = \mathbb{1}_{(1/8,3/8)}(x) + \mathbb{1}_{(5/8,7/8)}(x), $$
where $\mathbb{1}_E$ denotes the characteristic function of the set $E \subset [0, 1]$, and where $u_3^0$ is deduced from $u_1^0$ and $u_2^0$ by the relation $u_3^0 = 1 - u_1^0 - u_2^0$. The time step is chosen to be constant and equal to $\Delta t = 10^{-5}$ and final time as $T = 0.5$. The spatial mesh is chosen to be a uniform grid of the interval $(0, 1)$ containing $N$ subintervals.

The values of the cross-diffusion coefficients are chosen to be

$$c_{12} = c_{21} = 0.2, c_{13} = c_{31} = 1.0, c_{23} = c_{32} = 0.1, c^* = 0.1.$$

Figure 1 illustrates the evolution of the $L^1$ in time and space error of the approximate discrete solution as a function of $N$ (which is computed in comparison with an approximate solution computed on a very fine grid with $N_{\text{ref}} = 10^4$ cells).

<figure><img src="image_3.png" /><figcaption>FIGURE 1. Evolution of the $L^1$ space time error of the approximate solution as a function of the spatial discretization parameter.</figcaption></figure>

We numerically observe that the error decays like $O(1/N^2)$, in other words, showing that the scheme is second order accurate in space.

## 5.2. Two-dimensional test case.
We present here a two-dimensional test case. The number of species is kept to be $n = 3$ and the values of the cross-diffusion coefficients are now given by

$$(69) \qquad c_{12} = c_{21} = 0.1, c_{13} = c_{31} = 0.2, c_{23} = c_{32} = 2, c^* = 0.1.$$

The spatial domain $\Omega = (0, 1)^2$ is discretized using a cartesian uniform grid containing 70 cells in each direction. Time step is chosen to be $\Delta t = 10^{-5}$.

Figure 2 (respectively Figure 3 and Figure 4) shows the values of the concentration profiles $u_1, u_2, u_3$ at time $t = 0$ (respectively $t_1 = 8.5 \ 10^{-5}$ and $t_2 = 1 \ 10^{-3}$). Since the initial concentration profiles vanish on some parts of the domain, the projection step (66) plays an important role. A bare Newton-Raphson algorithm would produce undershoots of magnitude 0.03 for the first time step $\Delta t_1 = 5.10^{-6}$.

Since the coefficients $c_{12}$ and $c_{13}$ are much smaller than $c_{23}$, the initial interfaces between the different species are easily diffused for early times. Recall that $c_{ij}$ is an inverse diffusion coefficient. On Figure 4, one clearly sees that species 2 and 3 have difficulties to interdiffuse due to the high value of $c_{23}$, so that the species 2 remains essentially confined in a region where $u_3$ is small.

Our last figure is there to highlight both the decay of the discrete entropy and the exponential convergence towards equilibrium of the approximate solution. The exponential convergence in the continuous case was established in [41] thanks to a
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Initial profiles of the volume fractions.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 3. Profiles of the volume fractions at $t_1 = 8.5 \ 10^{-5}$.</figcaption></figure>

Logarithmic Sobolev inequality. A discrete counterpart of this inequality has been proved in [9], allowing to show the exponential convergence of the approximate solution towards the constant in space equilibrium following the lines of [41]. We omit the proof here and rather provide a numerical evidence.

Define $u^\infty := (u_{i,K}^\infty)_{i,K} \in (\mathbb{R}^T)^n$ by

$$u_{i,K}^\infty = \frac{1}{m_\Omega} M_i, \quad \forall 1 \le i \le n, \quad \forall K \in \mathcal{T},$$
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Profiles of the volume fractions at $t_2 = 1 \cdot 10^{-3}$.</figcaption></figure>

and $M_i$ is defined by (9), and by

$$H_{\mathcal{T}}(u^p|u^{\infty}) = \sum_{K \in \mathcal{T}} \sum_{i=1}^{n} m_K u_{i,K}^p \log \left( \frac{u_{i,K}^p}{u_{i,K}^{\infty}} \right) = E_{\mathcal{T}}(u^p) - E_{\mathcal{T}}(u^{\infty}) \geq 0$$

the relative entropy between the approximate solution $u^p$ at the $p^{th}$ time step and the long-time limit of $u$. Figure 5 shows that our approximate solution converges exponentially fast towards the correct long-time limit. The exponential convergence in $L^1$ can then be deduced from a Csiszár-Kullback inequality.

<figure><img src="image_7.png" /><figcaption>FIGURE 5. Evolution of the relative entropy $H_{\mathcal{T}}(u^p|u^{\infty})$ as a function of time.</figcaption></figure>

**Acknowledgements.** The authors acknowledge support from project COMODO (ANR-19-CE46-0002). CC also acknowledges support from Labex CEMPI (ANR-11-LABX-0007-01).
REFERENCES

[1] A. Ait Hammou Oulhaj. Numerical analysis of a finite volume scheme for a seawater intrusion model with cross-diffusion in an unconfined aquifer. *Numer. Methods Partial Differential Equations*, 34(3):857–880, 2018.

[2] L. Almeida, F. Bubba, B. Perthame, and C. Pouchol. Energy and implicit discretization of the Fokker-Planck and Keller-Segel type equations. *Netw. Heterog. Media*, 14(1):23–41, 2019.

[3] B. Andreianov, M. Bendahmane, and R. Ruiz-Baier. Analysis of a finite volume method for a cross-diffusion model in population dynamics. *Math. Models Methods Appl. Sci.*, 21(2):307–344, 2011.

[4] B. Andreianov, C. Cancès, and A. Moussa. A nonlinear time compactness result and applications to discretization of degenerate parabolic–elliptic PDEs. *J. Funct. Anal.*, 273(12):3633–3670, 2017.

[5] F. R. A. Aznaran, P. E. Farrell, C. W. Monroe, and A. J. Van-Brunt. Finite element methods for multicomponent convection-diffusion. arXiv:2208.11949, 2022.

[6] J.W. Barrett and J.F. Blowey. Finite element approximation of a nonlinear cross-diffusion population model. *Numer. Math.*, 98(2):195–221, 2004.

[7] M. Bessemoulin-Chatard. A finite volume scheme for convection-diffusion equations with nonlinear diffusion derived from the Scharfetter-Gummel scheme. *Numer. Math.*, 121(4):637–670, 2012.

[8] M. Bessemoulin-Chatard, C. Chainais-Hillairet, and M.-H. Vignal. Study of a finite volume scheme for the drift-diffusion system. Asymptotic behavior in the quasi-neutral limit. *SIAM J. Numer. Anal.*, 52(4):1666–1691, 2014.

[9] M. Bessemoulin-Chatard and A. Jüngel. A finite volume scheme for a Keller-Segel model with additional cross-diffusion. *IMA J. Numer. Anal.*, 34(1):96–122, 2014.

[10] D. Bothe. On the Maxwell-Stefan approach to multicomponent diffusion. In J. Escher, P. Guidotti, M. Hieber, P. Mucha, J. W. Prüss, Y. Shibata, G. Simonett, C. Walker, and W. Zajaczkowski, editors, *Parabolic Problems: The Herbert Amann Festschrift*, pages 81–93. Springer Basel, Basel, 2011.

[11] L. Boudin, B. Grec, and F. Salvarani. A mathematical and numerical analysis of the Maxwell-Stefan diffusion equations. *Discr. Cont. Dyn. Syst. B*, 17(5):427–440, 2012.

[12] L. Boudin, B. Grec, and F. Salvarani. The Maxwell-Stefan diffusion limit for a kinetic model of mixtures. *Acta Appl. Math.*, 136:79–90, 2015.

[13] M. Braukhoff, I. Perugia, and P. Stocker. An entropy structure preserving space-time formulation for cross-diffusion systems: analysis and Galerkin discretization. *SIAM J. Numer. Anal.*, 60:364–395, 2022.

[14] M. Braukhoff, C. Raithel, and N. Zamponi. Partial Hölder regularity for solutions of a class of cross-diffusion systems with entropy structure. arXiv:2204.06080, 2022.

[15] C. Cancès. Energy stable numerical methods for porous media flow type problems. *Oil & Gas Science and Technology-Rev. IFPEN*, 73:1–18, 2018.

[16] C. Cancès, C. Chainais-Hillairet, J. Fuhrmann, and B. Gaudeul. A numerical analysis focused comparison of several finite volume schemes for a unipolar degenerated drift-diffusion model. HAL: hal-02194604, 2019.

[17] C. Cancès, C. Chainais-Hillairet, A. Gerstenmayer, and A. Jüngel. Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport. *Numer. Methods Partial Differential Equations*, 35(2):545–575, 2019.

[18] C. Cancès and B. Gaudeul. A convergent entropy diminishing finite volume scheme for a cross-diffusion system. *SIAM J. Numer. Anal.*, 58(5):2684–2710, 2020.

[19] C. Cancès and C. Guichard. Convergence of a nonlinear entropy diminishing Control Volume Finite Element scheme for solving anisotropic degenerate parabolic equations. *Math. Comp.*, 85(298):549–580, 2016.

[20] C. Cancès and C. Guichard. Numerical analysis of a robust free energy diminishing finite volume scheme for parabolic equations with gradient structure. *Found. Comput. Math.*, 17(6):1525–1584, 2017.

[21] José A Carrillo, Francis Filbet, and Markus Schmidtchen. Convergence of a finite volume scheme for a system of interacting species with cross-diffusion. *Numerische Mathematik*, 145(3):473–511, 2020.
[22] C. Chainais-Hillairet. Entropy method and asymptotic behaviours of finite volume schemes. In *Finite volumes for complex applications. VII. Methods and theoretical aspects*, volume 77 of Springer Proc. Math. Stat., pages 17–35. Springer, Cham, 2014.

[23] E. S. Daus, A. Jüngel, and A. Zurek. Convergence of a finite-volume scheme for a degenerate-singular cross-diffusion system for biofilms. *IMA J. Numer. Anal.*, 41(2):935–973, 2021.

[24] K. Deimling. *Nonlinear functional analysis*. Courier Corporation, 2010.

[25] J. Droniou. Finite volume schemes for diffusion equations: introduction to and review of modern methods. *Math. Models Methods Appl. Sci.*, 24(08):1575–1619, 2014.

[26] V. Ehrlacher, C. Cancès, and L. Monasse. Finite volume scheme for the Stefan-Maxwell model. https://doi.org/10.5281/zenodo.3934286, July 2020.

[27] R. Eymard, T. Gallouët, C. Guichard, R. Herbin, and R. Masson. TP or not TP, that is the question. *Comput. Geosci.*, 18(3-4):285–296, 2014.

[28] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. Ciarlet, P. G. (ed.) et al., in Handbook of numerical analysis. North-Holland, Amsterdam, pp. 713–1020, 2000.

[29] M. Frittelli, A. Madzvamuse, I. Sgura, and C. Venkataraman. Lumped finite elements for reaction-cross-diffusion systems on stationary surfaces. *Computers & Mathematics with Applications*, 74(12):3008–3023, 2017.

[30] K. Gärtner and L. Kamenski. Why do we need Voronoi cells and Delaunay meshes? In *Numerical Geometry, Grid Generation and Scientific Computing*, pages 45–60. Springer, 2019.

[31] B. Gaudeul and J. Fuhrmann. Entropy and convergence analysis for two finite volume schemes for a Nernst-Planck-Poisson system with ion volume constraints. *Numer. Math.*, 151(1):99–149, 2022.

[32] A. Gerstenmayer and A. Jüngel. Comparison of a finite-element and finite-volume scheme for a degenerate cross-diffusion system for ion transport. *Comput. Appl. Math.*, 38(3):Art. 108, 23, 2019.

[33] V. Giovangigli. Multicomponent flow modeling. *Science China Mathematics*, 55(2):285–308, 2012.

[34] A. Gurusamy and K. Balachandran. Finite element method for solving Keller-Segel chemotaxis system with cross-diffusion. *International Journal of Dynamics and Control*, 6(2):539–549, 2018.

[35] R. Herbin. An error estimate for a finite volume scheme for a diffusion-convection problem on a triangular mesh. *Numer. Methods Partial Differential Equations*, 11(2):165–173, 1995.

[36] H. Huo, X. and Liu, A. E. Tzavaras, and S. Wang. An energy stable and positivity-preserving scheme for the Maxwell-Stefan diffusion system. *SIAM J. Numer. Anal.*, 59(5):2321–2345, 2021.

[37] X. Huo, A. Jüngel, and A. E. Tzavaras. Weak-strong uniqueness for Maxwell-Stefan systems. *SIAM J. Math. Anal.*, 54(3):3215–3252, 2022.

[38] A. Jüngel. The boundedness-by-entropy method for cross-diffusion systems. *Nonlinearity*, 28(6):1963, 2015.

[39] A. Jüngel. *Entropy methods for diffusive partial differential equations*. SpringerBriefs in Mathematics. Springer, [Cham], 2016.

[40] A. Jüngel and O. Leingang. Convergence of an implicit Euler Galerkin scheme for Poisson-Maxwell-Stefan systems. *Adv. Comput. Math.*, 45(3):1469–1498, 2019.

[41] A. Jüngel and I. V. Stelzer. Existence analysis of Maxwell-Stefan systems for multicomponent mixtures. *SIAM J. Math. Anal.*, 45(4):2421–2440, 2013.

[42] A. Jüngel and A. Zurek. A convergent structure-preserving finite-volume scheme for the Shigesada-Kawasaki-Teramoto population system. *SIAM J. Numer. Anal.*, 59(4):2286–2309, 2021.

[43] J. Kou, S. Sun, and X. Wang. Linearly decoupled energy-stable numerical methods for multicomponent two-phase compressible flow. *SIAM J. Numer. Anal.*, 56(6):3219–3248, 2018.

[44] J. Leray and J. Schauder. Topologie et équations fonctionnelles. *Ann. Sci. École Norm. Sup.*, 51((3)):45–78, 1934.

[45] J. C. Maxwell. IV. on the dynamical theory of gases. *Phil. Trans. Royal Soc. London*, 157:49–88, 1867.

[46] H. Murakawa. A linear finite volume method for nonlinear cross-diffusion systems. *Numer. Math.*, 136(1):1–26, 2017.
[47] K. S. C. Peerenboom, J. van Dijk, J. H. M. ten Thije Boonkkamp, L. Liu, W. J. Goedheer, and J. J. A. M. van der Mullen. Mass conservative finite volume discretization of the continuity equations in multi-component mixtures. *J. Comput. Phys.*, 230:3525–3537, 2011.

[48] J. Shen, J. Xu, and J. Yang. A new class of efficient and robust energy stable schemes for gradient flows. *SIAM Rev.*, 61(3):474–506, 2019.

[49] J. Stefan. Ueber das Gleichgewicht und Bewegung, insbesondere die Diffusion von Gemischen. *Sitzungsberichte der Kaiserlichen Akademie der Wissenschaften Wien, 2te Abteilung a*, 63, 1871.

[50] J. M. Stockie, K. Promislow, and B. R. Wetton. A finite volume method for multicomponent gas transport in a porous fuel cell electrode. *Int. J. Numer. Meth. Fluids*, 41:577–599, 2003.

[51] Z. Sun, J. A. Carrillo, and C.-W. Shu. A discontinuous Galerkin method for nonlinear parabolic equations and gradient flow problems with interaction potentials. *J. Comput. Phys.*, 352:76–104, 2018.

[52] Z. Sun, J. A. Carrillo, and C.-W. Shu. An entropy stable high-order discontinuous Galerkin method for cross-diffusion gradient flow systems. *Kinet. Relat. Models*, 12(4):885–908, 2019.

[53] A. Van-Brunt, P. E. Farrell, and C. W. Monroe. Augmented saddle-point formulation of the steady-state Stefan-Maxwell diffusion problem. *IMA J. Numer. Anal.*, 10 2021. drab067.

[54] J. A. Wesselingh, Rajamani Krishna, et al. *Mass transfer in multicomponent mixtures*. Delft University Press Delft, 2000.

CLÉMENT CANCÈS (clement.cances@inria.fr), INRIA, UNIV. LILLE, CNRS, UMR 8524 - LABORATOIRE PAUL PAINLEVÉ, F-59000 LILLE.

VIRGINIE EHRLACHER (virginie.ehrlacher@enpc.fr), CERMICS, ECOLE DES PONTS PARIS-TECH AND INRIA PARIS, UNIVERSITÉ PARIS-EST, 6-8 AVENUE BLAISE PASCAL, 77455, MARNE-LA-VALLÉE, FRANCE.

LAURENT MONASSE (laurent.monasse@inria.fr), UNIVERSITÉ CÔTE D'AZUR, INRIA, CNRS, LABORATOIRE J.A. DIEUDONNÉ, TEAM COFFEE, PARC VALROSE, 06108 NICE CEDEX 02, FRANCE.