# Finite Volumes for a dissipative free boundary problem

Clément Cancès, Claire Chainais-Hillairet, Amélie Dupouy

► To cite this version:

Clément Cancès, Claire Chainais-Hillairet, Amélie Dupouy. Finite Volumes for a dissipative free boundary problem. 2025. hal-05266770

HAL Id: hal-05266770

https://hal.science/hal-05266770v1

Preprint submitted on 18 Sep 2025

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 A DISSIPATIVE FREE BOUNDARY PROBLEM

CLÉMENT CANCÈS, CLAIRE CHAINAIS-HILLAIRET, AND AMÉLIE DUPOUY

**ABSTRACT.** We study a toy model for the evolution of the oxygen concentration in an oxide layer. It consists in a transient convection diffusion equation in a one-dimensional domain of variable width. The motions of the boundaries are governed by the traces of the concentration. We exhibit a necessary and sufficient condition on the parameters involved in the model for the existence of a unique traveling-wave solution. Moreover, we show that the model admits some universal entropy structure, in the sense that any convex function of the concentration yields a dissipated free energy (up to exchanges with the outer environment at the boundaries). We propose then an implicit in time arbitrary Lagrangian-Eulerian finite volume scheme based on Scharfetter-Gummel fluxes. It is shown to be unconditionally convergent, to preserve exactly the travelling wave, and to dissipate all the aforementioned free energies. Numerical experiments show that our scheme is first order accurate in time and second order in space, and that the transient solution converges in the long-time limit towards the traveling-wave solution.

## 1. INTRODUCTION

1.1. **Motivation and context.** The so-called *Diffusion Poisson Coupled Model* (DPCM) introduced in [4] aims at understanding finely the dynamics of an oxide (magnetite) layer at the surface of iron canisters to be stored in deep underground nuclear waste repositories. Because of its complexity –3 different charged species evolving in a free boundary domain with retroaction of a self-consistent electric potential– and of its lack of mathematical structure –no Lyapunov functional could be identified so-far–, no satisfying well-posedness was obtained so far on the full model. Moreover, numerical experiments presented in [4, 5] show that, at least for particular choices of the parameters, the transient system converges in the long-time limit towards a travelling wave profile.

Global in time existence results in the simpler case where the boundaries are fixed are provided in [15, 10], while the existence of a travelling wave profile is shown in [14] assuming the electroneutrality of the oxide. The existence of a travelling wave accounting for the coupling with a self-consistent electric field is also rigorously established in [7] thanks to computer assisted proof technics, but only for some particular set of parameters.

In this paper we study a simplified model with only one chemical species and where electrical interactions are neglected. The goal there is to take advantage of the apparent simplicity of the model to answer some questions that were out of reach in more complex models. Is the existence of a travelling wave profile automatic, or are there some required conditions on the parameters to ensure the existence of such a particular solution? Can we extract a Lyapunov functional providing enough control on the problem so that we can rigorously prove the existence of a solution to the system? Does the travelling wave profile act as an attractor in the long-time limit? Can one design a numerical scheme that preserves at the discrete level the key features of the continuous problem? And can one rigorously establish the convergence of the scheme?

<footnote n="1991">1991 Mathematics Subject Classification. 35R35, 65M08, 65M12.</footnote>
We show in what follows that, even in our simplified setting, conditions are required on the parameters to make the existence of a travelling wave profile possible. In our greatly simplified context, we even get a necessary and sufficient condition for the existence of the travelling wave profile, for which we derive an explicit expression.

Besides, the problem under study appears to have an infinite number of free energies, which all dissipate up to some exchange terms with the outer environment. Any evolution of the solution profile or of the interface location generate an entropy production, and even more: the evolution is driven by the dissipation of the (infinite family of) free energies. In mathematical words, the system we consider here is a generalized gradient flow, in the sense of [30], and more precisely a port-gradient system [31] since there are exchanges between the oxide and its surrounding environment. From a thermodynamical perspective, the fluxes are driven by the variations of the so-called chemical potentials as in [10], whereas the motions of the interfaces are governed by the pressure (or equivalently by the grand potential), in good accordance with the literature [32, 11, 29]. Excepted for a very particular choice of the parameters, the system cannot reach an equilibrium, in opposition to the closely related configurations studied in [33] or [9]. As the numerical evidences show, the system rather converges towards a travelling wave profile. Such travelling waves for convection diffusion problems are also of great interest in other contexts, like for instance in models for cell motility [34, 1].

In order to approximate the solutions to the continuous problem, we introduce an implicit in time arbitrary Lagrangian Eulerian (ALE) finite volume scheme (see for instance [19]). In our one-dimensional context, this approach is equivalent to the one proposed in [5, 16], where a change of variable is introduced to formulate the problem on the fixed spatial domain (0, 1). Even after this change of variable, the continuous problem remains of convection diffusion type. We employ Scharfetter-Gummel type fluxes [35, 27, 17]. These fluxes exactly capture the travelling wave, and dissipate all the convex free energies of the system. Building on these strong properties, we establish the convergence of the scheme towards a weak solution thanks to compactness arguments.

We highlight in numerical experiments that the scheme is second order accurate in space and first order in time as expected. Furthermore, we observe that, in the long-time asymptotic regime, the solution to the problem converges exponentially fast towards the travelling wave profile (provided that the latter exists). Unfortunately, we have been unable to rigorously prove this result so far. Indeed, the port-gradient nature of the problem makes usual approaches such as [22, 26, 3, 24, 25, 6] fail in our case. This will be the purpose of future work.

## 1.2. Presentation of the continuous model.
In the simplified model studied in this paper, one chemical species –that should be thought as oxygen in the context of corrosion– evolves inside the one-dimensional oxide layer represented by the interval $(X_0(t), X_1(t))$, see Figure 1. The oxide is surrounded by two different outer environments, modeling an aqueous solution for $x < X_0(t)$ and a metal subject to oxidation for $x > X_1(t)$. Oxygen, the concentration of which being denoted by $u$ in the oxide, can enter the oxide layer from the solution, diffuse within the oxide, and then transform the metal into oxide thanks to an oxidation process.

From a mathematical perspective, the problem amounts to find $(u, X_0, X_1)$, where $X_\ell : [0, T) \to \mathbb{R}$ for $\ell \in \{0, 1\}$, and

$$u: Q_T = \bigcup_{t \in [0, T)} (\{t\} \times [X_0(t), X_1(t)]) \to \mathbb{R}_{\ge 0}$$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Schematic representation of the oxide in its environment.</figcaption></figure>

with $T > 0$, such that, for any $(t, x)$ in $Q_T$, it holds

$$ (1a) \qquad \partial_t u(t, x) + \partial_x J(t, x) = 0, $$

$$ (1b) \qquad J(t, x) = -\partial_x u(t, x) + (1-R)X_1'(t)u(t, x), $$

$$ (1c) \qquad J(t, X_0(t)) - u(t, X_0(t))X_0'(t) = a - b u(t, X_0(t)), $$

$$ (1d) \qquad J(t, X_1(t)) - u(t, X_1(t))X_1'(t) = 0, $$

$$ (1e) \qquad X_0'(t) = \alpha_0 - \beta_0 u(t, X_0(t)) + (1-R)X_1'(t), $$

$$ (1f) \qquad X_1'(t) = -\alpha_1 + \beta_1 u(t, X_1(t)), $$

$$ (1g) \qquad u(0, x) = u^{\text{init}}(x), $$

$$ (1h) \qquad X_0(0) = 0, $$

$$ (1i) \qquad X_1(0) = L^0, $$

The first equation (1a) is the conservation law on the concentration of oxygen inside the oxide. The flux $J$ in (1b) is composed of two terms: a diffusion flux, and a transport term related to the movement of the oxide with respect to the reference frame attached to the metal. The frame of the oxide moves in this reference frame with velocity

$$ v_{\text{ox}}(t) = (1-R)X_1'(t), $$

where $R$ is the so-called Pilling-Bedworth ratio encoding the difference in molar volumes between the metal and oxide. $R$ is assumed to be a given positive constant.

The conservative fluxes $\mathcal{F}_\ell$ at the interfaces $X_\ell$ are defined by

$$ (2) \qquad \mathcal{F}_\ell(t) = J(t, X_\ell(t)) - u(t, X_\ell(t))X_\ell'(t) \quad \text{for } \ell \in \{0, 1\}. $$

There is no oxygen exchange between the oxide and the metal, which transmits through the condition on the conservative fluxes in (1d). The other condition (1c) is obtained thanks to the mass action law, that becomes in this case a Robin condition at the interface between the solution and the oxide. In (1c), the kinetic parameters $a$ and $b$ are given and positive.

Both (1e) and (1f) describe the evolution of the oxide boundaries along time. They involve the positive given kinetic parameters $\alpha_0, \beta_0, \alpha_1$, and $\beta_1$. The growth of the oxide over the metal is given by (1f), which comes from the mass action law, and (1e) combines both the oxide frame velocity $v_{\text{ox}}$, and the dissolution velocity, denoted $v_d$, i.e.

$$ X_0'(t) = v_{\text{ox}}(t) + v_d(t), \quad \text{with} \quad v_d(t) = \alpha_0 - \beta_0 u(t, X_0(t)). $$

Throughout the paper, we denote by $L(t) = X_1(t) - X_0(t) > 0$ the thickness of the oxide. Since the problem is translation invariant, the initial value of the solution-oxide layer interface is set to 0 for simplicity in (1h). The initial position of the interface $X_1$ prescribed by (1i) then directly
encodes the initial thickness $L^0 > 0$ of the oxide layer. Finally, the initial concentration profile $u^{\text{init}}$ appearing in (1g) is assumed to be positive and bounded.

1.3. **Presentation of the numerical scheme.** Before defining our numerical scheme, let us first introduce some notation. We set a finite time horizon $T > 0$, the value of which will be discussed later on. We also introduce a mesh of $[0, 1]$, composed of $I > 0$ cells $(\xi_{i-\frac{1}{2}}, \xi_{i+\frac{1}{2}})$ such that

$$0 = \xi_{\frac{1}{2}} < \xi_{\frac{3}{2}} < \dots < \xi_{I-\frac{1}{2}} < \xi_{I+\frac{1}{2}} = 1.$$

The size of each cell is denoted $h_i = \xi_{i+\frac{1}{2}} - \xi_{i-\frac{1}{2}}$ for $1 \le i \le I$, and the cell centers $\xi_i$, $0 \le i \le I$, are defined as

$$\xi_i = \frac{\xi_{i-\frac{1}{2}} + \xi_{i+\frac{1}{2}}}{2} \text{ for } 1 \le i \le I, \quad \xi_0 = \xi_{\frac{1}{2}}, \text{ and } \xi_{I+1} = \xi_{I+\frac{1}{2}}.$$

The distance between two consecutive cell centers $\xi_i$ is denoted by $h_{i+\frac{1}{2}} = \xi_{i+1} - \xi_i$ for $0 \le i \le I+1$. The mesh size $h$ is then defined as $h = \max_{1 \le i \le I} h_i$.

For notation convenience, we consider uniform time discretization. The time step is denoted $\Delta t$, and it is assumed that there exists an integer $N$ such that $N\Delta t = T$. The sequence $(t_n)_{0 \le n \le N}$ is defined as $t_n = n\Delta t$ for each $0 \le n \le N$.

Based on the (0, 1) mesh, at each set time $t > 0$, a mesh of $(X_0(t), X_1(t))$ is defined, with $x_{i+\frac{1}{2}}(t) = X_0(t) + L(t)\xi_{i+\frac{1}{2}}$.

It is then possible to look at the evolution of the oxygen concentration in a mesh cell $[x_{i-\frac{1}{2}}(t), x_{i+\frac{1}{2}}(t)]$, with $1 \le i \le I$. Using (1a), it comes that

$$(3) \qquad \frac{d}{dt} \int_{x_{i-\frac{1}{2}}(t)}^{x_{i+\frac{1}{2}}(t)} u(t, x)dx = -\mathcal{F}_{i+\frac{1}{2}}(t) + \mathcal{F}_{i-\frac{1}{2}}(t),$$

where $\mathcal{F}_{i+\frac{1}{2}}(t) = J(t, x_{i+\frac{1}{2}}(t)) - u(x_{i+\frac{1}{2}}(t))x'_{i+\frac{1}{2}}(t)$. Bearing in mind the expression (1b) of the flux $J$, we obtain

$$\mathcal{F}_{i+\frac{1}{2}}(t) = -\partial_x u(t, x_{i+\frac{1}{2}}(t)) + u(t, x_{i+\frac{1}{2}}(t))v_{i+\frac{1}{2}}(t),$$

with the velocity $v_{i+\frac{1}{2}}$ being defined by

$$(4) \qquad v_{i+\frac{1}{2}}(t) = (1-R)X_1'(t) - L'(t)\xi_{i+\frac{1}{2}} - X_0'(t).$$

The initial conditions (1g)-(1i) are discretized into

$$(5) \qquad X_0^0 = 0, \quad X_1^0 = L^0, \quad \text{and} \quad u_i^0 = \frac{1}{L^0 h_i} \int_{L^0 \xi_{i-\frac{1}{2}}}^{L^0 \xi_{i+\frac{1}{2}}} u^{\text{init}}(x)dx.$$
A finite volume method is used to approximate the left-hand side of (3), and an implicit Euler method is used for the time discretization. This leads to the following numerical scheme:

$$ (6a) \qquad h_i \frac{L^n u_i^n - L^{n-1} u_i^{n-1}}{\Delta t} + \mathcal{F}_{i+\frac{1}{2}}^n - \mathcal{F}_{i-\frac{1}{2}}^n = 0, \quad \text{for } 1 \le i \le I, $$

$$ (6b) \qquad \mathcal{F}_{\frac{1}{2}}^n - a + bu_0^n = 0, $$

$$ (6c) \qquad \mathcal{F}_{I+\frac{1}{2}}^n = 0, $$

$$ (6d) \qquad \frac{X_0^n - X_0^{n-1}}{\Delta t} - \alpha_0 + \beta_0 u_0^n - (1-R) \frac{X_1^n - X_1^{n-1}}{\Delta t} = 0, $$

$$ (6e) \qquad \frac{X_1^n - X_1^{n-1}}{\Delta t} + \alpha_1 - \beta_1 u_{I+1}^n = 0, $$

$$ (6f) \qquad L^n - X_1^n + X_0^n = 0. $$

As already mentioned in the introduction, we use Scharfetter-Gummel fluxes [35], that involve the so-called Bernoulli function defined by $B(r) = \frac{r}{e^r - 1}$, with continuous extension $B(0) = 1$, in their definition:

$$ (7) \qquad \mathcal{F}_{i+\frac{1}{2}}^n = \frac{1}{L^n h_{i+\frac{1}{2}}} \left( B \left( -L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n \right) u_i^n - B \left( L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n \right) u_{i+1}^n \right). $$

In the above formula, the discrete velocities $v_{i+\frac{1}{2}}^n$ are natural discrete counterparts of the semi-discrete quantities (4). Introducing the notation $\delta^n[f] = \frac{f^n - f^{n-1}}{\Delta t}$ for $f = X_0, X_1, L$, they write

$$ (8) \qquad v_{i+\frac{1}{2}}^n = (1-R)\delta^n[X_1] - \xi_{i+\frac{1}{2}}\delta^n[L] - \delta^n[X_0]. $$

The choice of the specific expression (7) for the numerical fluxes is motivated mainly by the fact that it guarantees the preservation of steady-state profiles such as the travelling wave in the case of this paper (see Section 3.1), while reaching the optimal second order accuracy in space [28]. More specifically, the Bernoulli function has some computational properties listed below that will be useful in the rest of the paper.

**Lemma 1.1.** The Bernoulli function is 1-Lipschitz continuous, positive, decreasing and satisfies:

(i) For $r \in \mathbb{R}$, $B(-r) - B(r) = r$;

(ii) For $r \in \mathbb{R}$, $B(r)e^r = B(-r)$;

(iii) For $\theta \in [0, 1]$, $p, q, r \in \mathbb{R}$:

$$ (9) \qquad B(-r)p - B(r)q = (\theta B(r) + (1-\theta)B(-r))(p-q) + r((1-\theta)q + \theta p). $$

1.4. **Outline of the paper.** The aim of this paper is to conduct both a theoretical and a numerical analysis of the simplified corrosion model (1). One main goal is to prove the existence of a solution to the system (1), which is obtained thanks to the numerical analysis of the scheme (6).

In a first part, the existence and the uniqueness of a travelling wave profile, as presented in [14], is established under some necessary and sufficient condition on the parameters of the model. Then, a result of decreasing free energy is established in Proposition 2.4, and leads to several bounds on any solution of (1) (see Proposition 2.5).

Then, the analysis of the numerical scheme (6) is conducted, beginning with the establishment of Proposition 3.2, which is a discrete counterpart of Proposition 2.4. One can then prove the same bounds as for the continuous case with an analogous reasoning. This, combined with bounds on the domain size (see Proposition 3.5) leads to the existence of a solution to the numerical scheme
(1) (Theorem 3.7). Some numerical experiments are provided in Section 5 to illustrate the behavior of an approximate solution on several test cases. The numerical simulations also illustrate the long-time behaviour of the approximate solution, that numerically converges to a travelling wave profile, recovering the results stated in [4] and [5].

The main result of the paper is the convergence of the solution of the scheme (6) towards a weak solution of the system (1) (Theorem 4.7), and it is the goal of Section 4.

## 2. SOME PROPERTIES OF THE CONTINUOUS MODEL

In this section, the goal is to study the continuous model (1). First, the notion of travelling wave is defined, and its existence is proved in Proposition 2.1. In a second part, the non-decreasing behaviour of free energies associated to the problem (see Proposition 2.4) leads to various energy estimates, particularly $L^\infty$ bounds, that are presented in Proposition 2.5.

### 2.1. Travelling wave profile.
The existence of a particular solution, that can be seen as a non-equilibrium steady state (as in [31] for example), is presented in this subsection. The existence of this particular profile, called a travelling wave profile, has been studied in [7] for the general DPCM ([4]). In the case of the model (1), the notion of travelling wave is defined as follows.

**Definition 1.** A triplet $(u, X_0, X_1)$ is called a travelling wave solution to (1) if there exist a profile $\hat{u}: \mathbb{R} \to \mathbb{R}_{\ge 0}$, a velocity $\hat{c} \in \mathbb{R}$ and a width $\hat{L} \in \mathbb{R}_{>0}$ such that for all $t \in \mathbb{R}_{\ge 0}$

$$ (10) \qquad X_0'(t) = X_1'(t) = \hat{c} \quad \text{and} \quad X_1(t) - X_0(t) = \hat{L}, \\ \text{and for all } (t, x) \in \mathbb{R}_{\ge 0} \times \bigcup_{t \in \mathbb{R}_{\ge 0}} (X_0(t), X_1(t)) \\ (11) \qquad \qquad \qquad \qquad \qquad \qquad u(t, x) = \hat{u}(x - \hat{c}t). $$

**Proposition 2.1.** Assume that the set of parameters satisfies either

$$ (12) \qquad \frac{\alpha_0 + R\alpha_1}{\beta_0 + R\beta_1} < \frac{a}{b} < \frac{\alpha_0}{\beta_0}, \\ \text{or} \\ (13) \qquad \frac{\alpha_0 + R\alpha_1}{\beta_0 + R\beta_1} > \frac{a}{b} > \frac{\alpha_0}{\beta_0}. $$

Then, there exists a unique travelling wave solution to (1), explicitly given by

$$ (14a) \qquad \hat{u}(y) = \frac{a}{b}e^{-R\hat{c}y} \quad \text{for } y \in [0, \hat{L}], \\ (14b) \qquad \hat{c} = \frac{1}{R}\left(\alpha_0 - \beta_0 \frac{a}{b}\right), \\ (14c) \qquad \hat{L} = -\frac{1}{R\hat{c}}\log\left(\frac{b}{\beta_1 a}(\alpha_1 + \hat{c})\right). $$

Under assumption

$$ (15) \qquad \qquad \qquad \qquad \qquad \qquad \qquad \frac{a}{b} = \frac{\alpha_0}{\beta_0} = \frac{\alpha_1}{\beta_1}, $$

there exists an infinite number of constant steady solutions with $\hat{c} = 0$, $\hat{u}(y) = \frac{a}{b}$ for all $y \in [0, \hat{L}]$ with $\hat{L} > 0$ being arbitrary. If neither (12), (13) nor (15) hold true, then there exists no travelling wave solution.
Proof. Suppose that there exists a solution $(u, X_0, X_1)$ that satisfies (10)-(11). The equation (1a) can then be rewritten in terms of $\hat{u}$ as

$$-\partial_y^2 \hat{u} - R\hat{c}\partial_y \hat{u} = \partial_y \hat{\mathcal{F}} = 0,$$

where we have set $\hat{\mathcal{F}} = -\partial_y \hat{u} - R\hat{c}\hat{u}$. Since the boundary condition (1d) at $X_1$ gives $\hat{\mathcal{F}}(X_1(t)) = 0$ for any $t$, it leads to $\hat{\mathcal{F}} = 0$. We further deduce from the other boundary condition (1c) that $\hat{u}(0) = \frac{a}{b}$. Hence, the travelling wave profile is given by (14a).

Plugging this expression in (1e) gives exactly the relation (14b) for $\hat{c}$. We note that the sign of $\hat{c}$ coincides with the sign of $\frac{\alpha_0}{\beta_0} - \frac{a}{b}$.

Assume first that $\hat{c} = 0$, or equivalently that $\frac{a}{b} = \frac{\alpha_0}{\beta_0}$, then we deduce from (14a) that $\hat{u}(y) = \frac{a}{b}$ for all $y \in [0, \hat{L}]$. This is compatible with (1f) if and only if (15) holds. Under this latter condition, the size of the domain $\hat{L}$ can be chosen arbitrarily as long as it remains positive. In other words, there exists an infinite number of constant equilibria under condition (15).

Assume now that $\hat{c} \neq 0$. We infer from (1f) that $\frac{b(\alpha_1 + \hat{c})}{a\beta_1} = e^{-R\hat{c}\hat{L}}$. This relation can be inverted into (14c) if and only if the left-hand side is positive, i.e. if $\alpha_1 + \hat{c} > 0$. This property is automatically satisfied if $\hat{c} > 0$, i.e. when $\frac{\alpha_0}{\beta_0} > \frac{a}{b}$. But formula (14c) gives a positive thickness $\hat{L}$ if and only if

$$\frac{b(\alpha_1 + \hat{c})}{\beta_1 a} < 1 \Leftrightarrow \frac{\alpha_0 + R\alpha_1}{\beta_0 + R\beta_1} < \frac{a}{b}.$$

Similarly, in the case $\hat{c} < 0$, i.e. when $\frac{\alpha_0}{\beta_0} < \frac{a}{b}$, the thickness $\hat{L}$ provided by formula (14c) is positive if and only if

$$\frac{b(\alpha_1 + \hat{c})}{\beta_1 a} > 1 \Leftrightarrow \frac{\alpha_0 + R\alpha_1}{\beta_0 + R\beta_1} > \frac{a}{b}.$$

This last condition is in this case more stringent than $\alpha_1 + \hat{c} > 0$, so that we end up with the simple conditions (12) or (13) on the parameters.

To conclude this proof, one readily check that, when (12) or (13) hold true, then $u(t, x) = \hat{u}(x - \hat{c}t)$, $X_0(t) = \hat{c}t$ and $X_1(t) = \hat{L} + \hat{c}t$ with $\hat{u}, \hat{c}$ and $\hat{L}$ given by (14) is a solution to the continuous problem (1). $\square$

**Remark 2.2.** In the rest of the paper, the study is carried out under assumption (12), which corresponds to the physical configuration encountered in the corrosion process, where metal is transformed into oxide, which itself dissolves in the solution. Moreover, remark that the assumption (12) immediately implies that $\frac{\alpha_1}{\beta_1} < \frac{\alpha_0}{\beta_0}$.

2.2. **Energy estimates.** In order to obtain some a priori estimates on the function $u$ solution to (1), we study the evolution in time of some free energies. In our case, any convex function can generate a free energy for the system. But, as the system is not isolated, in the sense that there are exchanges of mass and of volume with the outer environments, the free energy is not decreasing along time. To circumvent this difficulty, we follow the path proposed in [10] as we introduce in the following definition what we refer to as the total free energy, encompassing contributions from the outer environments.
**Definition 2.** Let $\phi \in C^2(\mathbb{R})$ be an arbitrary convex function. The free energy $\mathcal{H}_\phi$ associated to $\phi$ corresponding to the state $(u(t), X_0(t), X_1(t))$ is defined by

$$ (16) \qquad \mathcal{H}_{\phi}(t) = \int_{X_0(t)}^{X_1(t)} \phi(u(t, x))dx, \quad t \ge 0. $$

Moreover, defining the pressure function $\pi_\phi$ by $\pi_\phi(r) = r\phi'(r) - \phi(r)$ for $r \in \mathbb{R}$, the total free energy associated to $\phi$ is defined as

$$ (17) \qquad \mathcal{H}_{\phi}^{\text{tot}}(t) = \mathcal{H}_{\phi}(t) + R\pi_{\phi}\left(\frac{\alpha_1}{\beta_1}\right)(X_1(t) - X_1(0)) \\ \qquad \qquad - \pi_{\phi}\left(\frac{\alpha_0}{\beta_0}\right)\int_{0}^{t}(\alpha_0 - \beta_0u(t, X_0(t))) - \phi'\left(\frac{a}{b}\right)\int_{0}^{t}(a - bu(t, X_0(t))). $$

**Remark 2.3.** Since $\pi'_\phi(r) = r\phi''(r)$ with $\phi$ convex, $\pi_\phi$ is non-decreasing over $\mathbb{R}_{\ge 0}$.

The following proposition is proved with formal computations, assuming that $u \ge 0$, $X_0$ and $X_1$ exist and are regular enough.

**Proposition 2.4.** Assume that there exists a solution $(u, X_0, X_1)$ to the system (1) such that $u$ is non-negative and $X_1 - X_0$ positive and let $\phi \in C^2(\mathbb{R})$. Then, there exists a non-negative function $\mathcal{D}_\phi$ such that, for any $t > 0$,

$$ \frac{d}{dt} \mathcal{H}_{\phi}^{\text{tot}}(t) = -\mathcal{D}_{\phi}(t) \le 0. $$

*Proof.* Let $t > 0$. Compute the derivative of $\mathcal{H}_\phi(t)$ and use the equation (1a) to get

$$ \frac{d}{dt} \mathcal{H}_{\phi}(t) = - \int_{X_0(t)}^{X_1(t)} \phi'(u) \partial_x J + X_1'(t) \phi(u(t, X_1(t))) - X_0'(t) \phi(u(t, X_0(t))). $$

Integrating by parts and using the relation between the fluxes $J$ and $\mathcal{F}$ at the boundaries (2) yields

$$ \frac{d}{dt} \mathcal{H}_{\phi}(t) = - \int_{X_0(t)}^{X_1(t)} |\partial_x u|^2 \phi''(u) - R X_1'(t) \pi_{\phi}(u(t, X_1(t))) \\ + (X_0'(t) - (1-R)X_1'(t)) \pi_{\phi}(u(t, X_0(t))) + \mathcal{F}_0 \phi'(u(t, X_0(t))). $$

Then, using the equations on the moving boundaries (1e)-(1f), one gets that

$$ \frac{d}{dt} \mathcal{H}_{\phi}^{\text{tot}}(t) = -\mathcal{D}_{\phi}(t), $$

where the dissipation term $\mathcal{D}_\phi(t)$ is

$$ \mathcal{D}_{\phi}(t) = \int_{X_0(t)}^{X_1(t)} |\partial_x u|^2 \phi''(u) - R(\alpha_1 - \beta_1 u(t, X_1(t))) \left( \pi_{\phi}(u(t, X_1(t))) - \pi_{\phi}\left(\frac{\alpha_1}{\beta_1}\right) \right) \\ - (\alpha_0 - \beta_0 u(t, X_0(t))) \left( \pi_{\phi}(u(t, X_0(t))) - \pi_{\phi}\left(\frac{\alpha_0}{\beta_0}\right) \right) \\ - (a - bu(t, X_0(t))) \left( \phi'(u(t, X_0(t))) - \phi'\left(\frac{a}{b}\right) \right). $$

Since $\pi_\phi$ and $\phi'$ are non-decreasing over $\mathbb{R}_{\ge 0}$ (see Remark 2.3), it easily comes that all the terms in $\mathcal{D}_\phi(t)$ are non-negative, and so is $\mathcal{D}_\phi$. $\square$
This result gives various a priori bounds on any solution to the problem. The following proposition presents the a priori $L^\infty$ bound obtained thanks to the energy estimates.

**Proposition 2.5.** Assume that there exists a solution $(u, X_0, X_1)$ to the system (1) such that $u$ is non-negative and $X_1 - X_0$ positive. Then, under the assumption (12), $u$ satisfies $m \le u \le M$, where $m$ and $M$ are defined by

$$ (18) \qquad m = \min \left\{ \text{ess-inf} \, u^{\text{init}}, \frac{\alpha_1}{\beta_1} \right\}, \quad M = \max \left\{ \text{ess-sup} \, u^{\text{init}}, \frac{\alpha_0}{\beta_0} \right\}. $$

*Proof.* Take $\phi \in C^2(\mathbb{R})$ convex such that $\phi = 0$ on $[m, M]$, and $\phi > 0$ everywhere else. Then, $\phi'$ and $\pi_\phi$ are identically zero in $[m, M]$, thus Assumption (12) ensures that $\phi(r) = \phi'(r) = \pi_\phi(r) = 0$ for $r \in \left\{\frac{\alpha_1}{\beta_1}, \frac{a}{b}, \frac{\alpha_0}{\beta_0}\right\}$. Hence for all $t \in [0, T]$, there holds $\mathcal{H}_\phi^{\text{tot}}(t) = \mathcal{H}_\phi(t)$ and Proposition 2.4 gives that $\frac{d}{dt}\mathcal{H}_\phi \le 0$. Since $u^{\text{init}} \in [m, M]$, it comes that $\mathcal{H}_\phi(0) = 0$. Hence $\mathcal{H}_\phi(t) \le 0$ for $t \ge 0$. Since $\phi \ge 0$ the definition (16) of $\mathcal{H}_\phi$ ensures that $\mathcal{H}_\phi(t) \ge 0$, thus $\mathcal{H}_\phi(t) = 0$. Finally, $\phi$ being non-negative and continuous, it comes that $\phi(u) \equiv 0$, ensuring that $m \le u \le M$. $\square$

### 3. ANALYSIS ON A FIXED MESH

As a first contribution to the analysis of the numerical scheme, we establish some properties of the scheme on a given mesh, i.e. without letting the discretization parameters tend to 0.

**3.1. Preservation of the travelling wave profile.** In this part, it is proved that the travelling wave profile is a (quasi-)stationary state for the scheme.

**Proposition 3.1.** The numerical scheme (6) preserves the travelling wave profile in the following sense: for $0 \le i \le I+1$ and $0 \le n \le N$, denote $u_i^n = \hat{u}(\xi_i)$ with $\hat{u}$ defined in (11), set $X_1^n = \hat{L} + \hat{c}n\Delta t$ and $X_0^n = \hat{c}n\Delta t$. Then, the sequence $((u_i^n)_{0 \le i \le I+1}, X_0^n, X_1^n)_{0 \le n \le N}$ satisfies the scheme (6).

*Proof.* Let $0 \le i \le I+1$ and $1 \le n \le N$. Denote $u_i^n = \hat{u}(\xi_i)$, and set $X_1^n = \hat{L} + \hat{c}n\Delta t$ and $X_0^n = \hat{c}n\Delta t$. The definitions (14c), (14b) of $\hat{c}$ and $\hat{L}$ ensure that the equations on the boundaries in the scheme (6e) and (6d) are satisfied.

Remark that $v_{i+\frac{1}{2}}^n = -R\hat{c}$. Then, the property (ii) of the Bernoulli function given in Lemma 1.1, combined with the exponential profile $\hat{u}$, implies the cancellation of all the numerical fluxes. This yields (6a), (6b), (6c) and concludes the proof. $\square$

**3.2. Discrete energy estimates.** To begin with, one can prove a result of decreasing energy that can be seen as the discrete counterpart of Proposition 2.4. Thus, it is necessary to define a discrete free energy as follows.

**Definition 3.** Let $\phi \in C^2(\mathbb{R})$ an arbitrary convex function. Set $n \ge 0$. The discrete free energy associated to $\phi$ is given by

$$ \mathcal{H}_{\phi}^{n} = \sum_{i=1}^{I} L^{n} h_{i} \phi(u_{i}^{n}). $$
Moreover, the discrete total free energy associated to $\phi$ is

$$ (19) \qquad \mathcal{H}_{\phi}^{\text{tot},n} = \mathcal{H}_{\phi}^{n} + \pi_{\phi} \left(\frac{\alpha_0}{\beta_0}\right) \sum_{k=1}^{n} \Delta t(\alpha_0 - \beta_0 u_0^k) \\ \qquad \qquad + \phi' \left(\frac{a}{b}\right) \sum_{k=1}^{n} \Delta t(a - bu_0^k) + R\pi_{\phi} \left(\frac{\alpha_1}{\beta_1}\right) \sum_{k=1}^{n} \Delta t(\alpha_1 - \beta_1 u_{I+1}^k). $$

**Proposition 3.2.** Assume that the scheme (6) admits a solution $((u_i^n)_{0\le i\le I+1}, X_0^n, X_1^n)_{0\le n\le N}$ such that, for all $0 \le n \le N$, $L^n > 0$ and $u_i^n \ge 0$ for all $0 \le i \le I+1$. Then, for any convex function $\phi \in C^2(\mathbb{R})$, there exists a non-negative dissipation term $\mathcal{D}_{\phi}^n$ such that the discrete derivative of $\mathcal{H}_{\phi}^{\text{tot},n}$ satisfies

$$ (20) \qquad \frac{\mathcal{H}_{\phi}^{\text{tot},n} - \mathcal{H}_{\phi}^{\text{tot},n-1}}{\Delta t} + \mathcal{D}_{\phi}^{n} \le 0. $$

*Proof.* Let us introduce the application $A : (l, w) \mapsto l\phi\left(\frac{w}{l}\right)$ in order to rewrite

$$ \mathcal{H}_{\phi}^{n} = \sum_{i=1}^{I} h_{i} A(L^{n}, L^{n} u_{i}^{n}). $$

As $\phi$ is convex, $A$ is jointly convex and verifies $\nabla A(l, w) = (-\pi_{\phi}(\frac{w}{l}), \phi'(\frac{w}{l}))$. From an elementary convexity inequality, we deduce that

$$ (21) \qquad \frac{\mathcal{H}_{\phi}^{n} - \mathcal{H}_{\phi}^{n-1}}{\Delta t} \le Y^n + Z^n $$

with

$$ Y^n = \sum_{i=1}^{I} h_i \phi'(u_i^n) \frac{L^n u_i^n - L^{n-1} u_i^{n-1}}{\Delta t}, \quad Z^n = -\delta^n [L] \sum_{i=1}^{I} h_i \pi_{\phi}(u_i^n). $$

Using the numerical scheme (6a) in the expression $Y^n$, performing a discrete integration by parts and using the discrete boundary conditions (6b) and (6c) provides

$$ (22) \qquad Y^n = (a - bu_0^n)\phi'(u_0^n) + \sum_{i=0}^{I} (\phi'(u_{i+1}^n) - \phi'(u_i^n)) \mathcal{F}_{i+\frac{1}{2}}^n. $$

For any $0 \le i \le I$ and $n \ge 0$, the mean value formula and the link between $\phi$ and $\pi_{\phi}$ ensure the existence of $\theta_{i+\frac{1}{2}}^n \in [0, 1]$ such that

$$ (23) \qquad \pi_{\phi}(u_{i+1}^n) - \pi_{\phi}(u_i^n) = \left(\theta_{i+\frac{1}{2}}^n u_i^n + (1 - \theta_{i+\frac{1}{2}}^n) u_{i+1}^n\right) \left(\phi'(u_{i+1}^n) - \phi'(u_i^n)\right). $$

Then, using the third property of the Bernoulli function given in (9) in Lemma 1.1, it comes that

$$ (24) \qquad \sum_{i=0}^{I} (\phi'(u_{i+1}^n) - \phi'(u_i^n)) \mathcal{F}_{i+\frac{1}{2}}^n = -\mathcal{D}_{\phi}^{\text{bulk},n} + \sum_{i=0}^{I} v_{i+\frac{1}{2}}^n (\pi_{\phi}(u_{i+1}^n) - \pi_{\phi}(u_i^n)), $$
where the quantity $\mathcal{D}_{\phi}^{\text{bulk},n}$ is defined by

$$ (25) \quad \mathcal{D}_{\phi}^{\text{bulk},n} = \sum_{i=0}^{I} \left( B \left( L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n \right) \theta_{i+\frac{1}{2}}^n + B \left( -L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n \right) (1-\theta_{i+\frac{1}{2}}^n) \right) \times \frac{\phi'(u_{i+1}^n) - \phi'(u_i^n)}{L^n h_{i+\frac{1}{2}}} (u_{i+1}^n - u_i^n) \geq 0. $$

Thanks to a discrete integration by parts, the last term in (24) rewrites

$$ \sum_{i=0}^{I} v_{i+\frac{1}{2}}^{n} (\pi_{\phi}(u_{i+1}^{n}) - \pi_{\phi}(u_{i}^{n})) = -\sum_{i=1}^{I} (v_{i+\frac{1}{2}}^{n} - v_{i-\frac{1}{2}}^{n})\pi_{\phi}(u_{i}^{n}) + v_{I+\frac{1}{2}}^{n}\pi_{\phi}(u_{I+1}^{n}) - v_{\frac{1}{2}}^{n}\pi_{\phi}(u_{0}^{n}). $$

The definition (8) of the discrete velocities implies that $(v_{i+\frac{1}{2}}^n - v_{i-\frac{1}{2}}^n) = -\delta^n[L]h_i$ and also fix values of $v_{I+\frac{1}{2}}^n$ and $v_{I-\frac{1}{2}}^n$, so that

$$ (26) \quad \sum_{i=0}^{I} v_{i+\frac{1}{2}}^{n} (\pi_{\phi}(u_{i+1}^{n}) - \pi_{\phi}(u_{i}^{n})) = -Z^{n} - R\delta^{n}[X_1]\pi_{\phi}(u_{I+1}^{n}) - ((1-R)\delta^{n}[X_1] - \delta^{n}[X_0])\pi_{\phi}(u_{0}^{n}). $$

Then, using the equations on $\delta^n[X_0]$ and $\delta^n[X_1]$ in (6), we deduce from (21), (22), (24) and (26):

$$ \frac{\mathcal{H}_{\phi}^{n} - \mathcal{H}_{\phi}^{n-1}}{\Delta t} + \mathcal{D}_{\phi}^{\text{bulk},n} + (bu_{0}^{n} - a)\phi'(u_{0}^{n}) + (\beta_{0}u_{0}^{n} - \alpha_{0})\pi_{\phi}(u_{0}^{n}) + R(\beta_{1}u_{I+1}^{n} - \alpha_{1})\pi_{\phi}(u_{I+1}^{n}) \leq 0. $$

In order to recover the expression of $\mathcal{H}_{\phi}^{\text{tot},n}$ given in (19), it only remains to use the non-decreasing behavior of $\phi'$ and $\pi_{\phi}$ (see Remark 2.3). Then, the following term is non-negative:

$$ (27) \quad \mathcal{D}_{\phi}^{\text{bound},n} = (\beta_0 u_0^n - \alpha_0) \left( \pi_{\phi}(u_0^n) - \pi_{\phi}\left(\frac{\alpha_0}{\beta_0}\right) \right) + (bu_0^n - a) \left( \phi'(u_0^n) - \phi'\left(\frac{a}{b}\right) \right) + R(\beta_1 u_{I+1}^n - \alpha_1) \left( \pi_{\phi}(u_{I+1}^n) - \pi_{\phi}\left(\frac{\alpha_1}{\beta_1}\right) \right) \geq 0. $$

Finally, the inequality (20) holds with the dissipation term $\mathcal{D}_{\phi}^{n}$ defined as

$$ (28) \qquad \mathcal{D}_{\phi}^{n} = \mathcal{D}_{\phi}^{\text{bulk},n} + \mathcal{D}_{\phi}^{\text{bound},n}. $$

□

This result will be useful in the rest of the paper to get many a priori estimates on any solution to the numerical scheme, thanks to different choices of energy density function $\phi$.
### 3.2.1. Uniform bounds on the concentrations.
Let us establish the same a priori $L^\infty$ bounds as for the continuous problem.

**Proposition 3.3.** Assume that the scheme (6) admits a solution $((u_i^n)_{0\le i\le I+1}, X_0^n, X_1^n)_{0\le n\le N}$ such that, for all $0 \le n \le N$, $L^n > 0$ and $u_i^n \ge 0$ for all $0 \le i \le I+1$. Then, the solution to the scheme satisfies

$$ (29) \qquad m \le u_i^n \le M \quad \forall 0 \le i \le I+1, \forall 0 \le n \le N $$

where $m$ and $M$ are defined in (18).

To prove this result, we adapt to the discrete setting the proof of Proposition 2.5.

*Proof.* The result is shown by induction on the time step $n \ge 0$. The property (29) being clearly verified for $n=0$, set $n > 0$ and assume that for any $0 \le i \le I+1$, $u_i^{n-1}$ verifies (29).

As in the proof of Proposition 2.5, we take $\phi \in C^2(\mathbb{R})$ convex such that $\phi = 0$ on $[m, M]$, $\phi > 0$ and strictly convex everywhere else and we apply Proposition 3.2. For such a function $\phi$, one has $\mathcal{H}_\phi^{\text{tot},n} = \mathcal{H}_\phi^n$ and $\mathcal{H}_\phi^{\text{tot},n-1} = \mathcal{H}_\phi^{n-1}$. Moreover, $\mathcal{H}_\phi^{n-1} = 0$ thanks to the induction hypothesis. Then, one gets that

$$ \mathcal{H}_\phi^n + \Delta t \mathcal{D}_\phi^n \le 0. $$

Since both $\mathcal{H}_\phi^n$ and $\mathcal{D}_\phi^n$ are also both non-negative, we deduce that $\mathcal{H}_\phi^n = 0$, which yields $\phi(u_i^n) = 0$ and therefore $m \le u_i^n \le M$ for $1 \le i \le I$, and $\mathcal{D}_\phi^n = 0$, which yields the same bounds for $u_0^n$ and $u_{I+1}^n$. $\square$

One then readily deduces the following statement from (6d)-(6f) and (8).

**Corollary 3.4.** Under the assumptions of Proposition 3.3, there holds

$$ (30) \qquad -\alpha_1 + \beta_1 m \le \delta^n[X_1] \le -\alpha_1 + \beta_1 M, $$

$$ (31) \qquad -\alpha_0 - R\alpha_1 + m(\beta_0 + R\beta_1) \le \delta^n[L] \le -\alpha_0 - R\alpha_1 + M(\beta_0 + R\beta_1). $$

As a consequence, there exists some $v_\phi$ and $v_\phi^\sharp$ depending only on the data of the continuous problem such that

$$ (32) \qquad v_\phi \le v_{i+\frac{1}{2}}^n \le v_\phi^\sharp, \quad \forall i \in \{0, \dots, I\}. $$

### 3.2.2. Bounds on the domain size.
In this section, two different kinds of estimates are established on the domain size. The first one presented below in Proposition 3.5, is a non-uniform estimate with respect to the discretization parameters, but it is global in time. It will be useful to prove the existence of a solution to the scheme (6).

**Proposition 3.5.** Suppose that the hypothesis (12) is satisfied, and that there exists a solution $((u_i^n)_{0\le i\le I+1}, X_0^n, X_1^n)$ to the numerical scheme (6) with $L^n > 0$ and $u_i^n \ge 0$ for all $0 \le i \le I+1$. Then $L^n > \min(\frac{m}{M}L^{n-1}, \hat{L})$.

*Proof.* Let $n \ge 1$. Suppose for contradiction that $L^n \le \min(\frac{m}{M}L^{n-1}, \hat{L})$. This implies $\delta^n[L] < 0$ since Assumption (12) ensures that $m < M$. Furthermore, owing to the $L^\infty$-bounds (29), one has

$$ h_i \frac{L^n u_i^n - L^{n-1} u_i^{n-1}}{\Delta t} \le 0, \quad \forall 1 \le i \le I. $$

Using the scheme (6), it comes that

$$ (33) \qquad \mathcal{F}_{i-\frac{1}{2}}^n \le \mathcal{F}_{i+\frac{1}{2}}^n \le \mathcal{F}_{I+\frac{1}{2}}^n = 0, \quad \forall 1 \le i \le I. $$
In particular, $\mathcal{F}_{\frac{n}{2}} = a - bu_0^n \le 0$, that is $u_0^n \ge \frac{a}{b}$. This implies

$$ (34) \qquad v_{i+\frac{1}{2}}^n = -\alpha_0 + \beta_0 u_0^n - \delta^n[L]\xi_{i+\frac{1}{2}} \ge -\alpha_0 + \beta_0 \frac{a}{b} = -R\hat{c}. $$

Moreover, one infers from (33), from the definition (7) of the numerical fluxes, and from property (ii) of Lemma 1.1 that $u_{i+1}^n \ge u_i^n \exp(L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n)$, which yields

$$ u_{I+1}^n \ge \exp\left(\sum_{i=0}^{I} L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n\right) u_0^n \ge u_0^n e^{-L^n R \hat{c}} \ge \frac{a}{b} e^{-\hat{L} R \hat{c}} = \hat{u}(1). $$

The relation (6e) then provides that $\delta^n[X_1] \ge \hat{c}$. But, noticing that (6d) combined with (6f) rewrites

$$ \delta^n[L] = R\delta^n[X_1] - (\alpha_0 - \beta_0 u_0^n), $$

we then deduce that $\alpha_0 - \beta_0 u_0^n > R\hat{c}$, which contradicts $u_0^n \ge \frac{a}{b}$ obtained thanks to (33). Therefore, the starting hypothesis in this proof can not be satisfied, which concludes the proof of Proposition 3.5. $\square$

Now, the second bound presented in this section is uniform with respect to the discretization, but only up to some finite time horizon $T$. As it is a direct consequence of (31), we do not detail the proof here.

**Proposition 3.6.** Suppose that the hypothesis (12) is satisfied, and that there exists a solution $((u_i^n)_{0 \le i \le I+1}, X_0^n, X_1^n)$ to the numerical scheme (6) with $u_i^n \ge 0$ for all $0 \le i \le I+1$ for $T < \frac{L^0}{(\alpha_0 + R\alpha_1) - m(\beta_0 + R\beta_1)}$. Then, there exist two positive constants $\underline{L}$ and $\overline{L}$, depending only on $T, L^0, m, M$ and the kinetic parameters, such that for all $n \in \{1, \dots, N\}$, it holds

$$ (35) \qquad 0 < \underline{L} \le L^n \le \overline{L}. $$

### 3.3. Existence of a solution to the numerical scheme.
The main result of this section can finally be proved.

**Theorem 3.7.** Under the assumption (12), for a given mesh and a given time step, the numerical scheme (6) admits a solution $((u_i^n)_{0 \le i \le I+1}, X_0^n, X_1^n)$ for any $0 \le n \le N$.

*Proof.* The proof relies on a Brouwer topological degree argument (see for instance Chapter 1, Section 3.1 in [18]). Our goal is then to construct a homotopy to pass from an invertible linear system for $\lambda = 0$ to the system (6) for $\lambda = 1$ in such a way that the solution $((u_i^{(\lambda)})_{0 \le i \le I+1}, X_1^{(\lambda)}, L^{(\lambda)})$ for $\lambda \in [0, 1]$ remains in some open set $\Omega$ containing the solution for $\lambda = 0$ and such that the scheme is uniformly continuous on $\bar{\Omega}$. Note that we use $((u_i^{(\lambda)})_{0 \le i \le I+1}, X_1^{(\lambda)}, L^{(\lambda)})$ as main unknowns there, and that $X_0^{(\lambda)}$ is deduced a posteriori.
For $\lambda \in [0, 1]$, we seek $((u_i^{(\lambda)})_{0 \le i \le I+1}, X_1^{(\lambda)}, L^{(\lambda)})$ solution to

$$ 
\begin{aligned}
& \frac{\lambda(L^{(\lambda)} - L^{n-1})u_i^{(\lambda)} + L^{n-1}(u_i^{(\lambda)} - u_i^{n-1})}{\Delta t} h_i + \lambda (\mathcal{F}_{i+\frac{1}{2}}^{(\lambda)} - \mathcal{F}_{i-\frac{1}{2}}^{(\lambda)}) = 0, \quad 1 \le i \le I, \\
& \lambda (\mathcal{F}_{\frac{1}{2}}^{(\lambda)} - a + bu_0^{(\lambda)}) + (1-\lambda)(\beta_0 u_0^{(\lambda)} - \alpha_0) = 0, \\
& \lambda \mathcal{F}_{I+\frac{1}{2}}^{(\lambda)} + (1-\lambda)(\beta_1 u_{I+1}^{(\lambda)} - \alpha_1) = 0, \\
& R \frac{X_1^{(\lambda)} - X_1^{n-1}}{\Delta t} - \frac{L^{(\lambda)} - L^{n-1}}{\Delta t} - \alpha_0 + \beta_0 u_0^{(\lambda)} = 0, \\
& \frac{X_1^{(\lambda)} - X_1^{n-1}}{\Delta t} + \alpha_1 - \beta_1 u_{I+1}^{(\lambda)} = 0,
\end{aligned}
$$

with

$$ \mathcal{F}_{i+\frac{1}{2}}^{(\lambda)} = \frac{1}{L^{(\lambda)}h_{i+\frac{1}{2}}} \left( B(-L^{(\lambda)}h_{i+\frac{1}{2}}v_{i+\frac{1}{2}}^{(\lambda)})u_i^{(\lambda)} - B(L^{(\lambda)}h_{i+\frac{1}{2}}v_{i+\frac{1}{2}}^{(\lambda)})u_{i+1}^{(\lambda)} \right) $$

and

$$ v_{i+\frac{1}{2}}^{(\lambda)} = -R \frac{X_1^{(\lambda)} - X_1^{n-1}}{\Delta t} \xi_{i+\frac{1}{2}} + (1 - \xi_{i+\frac{1}{2}})(\beta_0 u_0^{(\lambda)} - \alpha_0). $$

For $\lambda = 0$, the unique solution to the above system is $u_i^{(0)} = u_i^{n-1}$ for $1 \le i \le I$, $u_0^{(0)} = \frac{\alpha_0}{\beta_0}$, $u_{I+1}^{(0)} = \frac{\alpha_1}{\beta_1}$, $X_1^{(0)} = X_1^{n-1}$ and $L^{(0)} = L^{n-1}$. Moreover, one can extend the proof of Proposition 3.2 and then Proposition 3.3, Corollary 3.4, and Proposition 3.5 to the case $\lambda \le 1$. In particular, one checks that

$$ 
\begin{aligned}
& m \le u_i^{(\lambda)} \le M, \\
& 0 < \min \left( \hat{L}, \frac{m}{M} L^{n-1} \right) \le L^{(\lambda)} \le L^{n-1} + \Delta t (M(\beta_0 + R\beta_1) - \alpha_0 - R\alpha_1), \\
& X_1^{n-1} + \Delta t (\beta_1 m - \alpha_1) \le X_1^{(\lambda)} \le X_1^{n-1} + \Delta t (\beta_1 M - \alpha_1).
\end{aligned}
$$

These uniform bounds allow to use Brouwer's topological degree to ensure that the system admits (at least) a solution for all $\lambda \in [0, 1]$. Note here that the positive lower bound on $L^{(\lambda)}$ is crucial to ensure the uniform continuity of the flux $\mathcal{F}_{i+\frac{1}{2}}^{(\lambda)}$. $\square$

### 4. CONVERGENCE OF THE SCHEME

The goal of this part is to establish rigorously, thanks to compactness and consistency arguments, that our scheme captures the solution to the continuous problem (1).

#### 4.1. Weak and approximate solutions.
As our goal is to establish the convergence of the discrete solution, the existence of which being guaranteed thanks to Theorem 3.7, one needs to first give a precise sense to the notion of solution to the continuous model (1). This is the purpose of the following definition.

**Definition 4.** A triplet $(u, X_0, X_1)$ is said to be a weak solution of (1) up to time $T > 0$ if it satisfies:

(i) $X_0$ and $X_1$ belong to $C^{0,1}([0, T], \mathbb{R})$, with $X_0(0) = 0$, $X_1(0) = L^0$ and $L(t) = X_1(t) - X_0(t) > 0$ for all $t \in [0, T)$;

(ii) $u \in L^\infty(Q_T)$, with $u \ge 0$ a.e. in $Q_T$, and $\partial_x u \in L^2(Q_T)$;
(iii) For any $\psi \in C_c^\infty(Q_T)$, there holds

$$ (36) \quad \int_{X_0(0)}^{X_1(0)} u^{\text{init}}\psi(0, \cdot) + \int_0^T \int_{X_0(t)}^{X_1(t)} (u\partial_t\psi + (-\partial_x u + (1-R)uX'_1)\partial_x\psi) \\ \qquad + \int_0^T (a - bu(\cdot, X_0(\cdot)))\psi(\cdot, X_0(\cdot)) = 0; $$

(iv) For a.e. $t \in [0, T)$, $X_0$ and $X_1$ satisfy

$$ X'_1(t) = -\alpha_1 + \beta_1 u_1(t) \quad \text{and} \quad X'_0(t) = \alpha_0 - \beta_0 u_0(t) + (1-R)X'_1(t), $$

where, for $\ell \in \{0, 1\}$, $u_\ell(t)$ denotes the trace of $u(t, \cdot)$ at $X_\ell(t)$.

**Remark 4.1.** Applying the change of variables

$$ (t, \xi) \in [0, T) \times [0, 1] \mapsto (t, x = X_0(t) + \xi L(t)) \in Q_T, $$

we obtain an equivalent definition of the weak solution. Indeed, $(u, X_0, X_1)$ is a weak solution of (1) in the sense of Definition 4 if and only if $X_0$ and $X_1$ satisfy (i) and there exists a positive function $w \in L^\infty([0, T) \times [0, 1]) \cap L^2([0, T] \times H^1(0, 1))$ such that

* For any $\varphi \in C_c^\infty([0, T) \times [0, 1])$, there holds

$$ (37) \quad \int_0^T \int_0^1 \left( L w \partial_t \varphi + \left( -\frac{1}{L} \partial_\xi w + ((1-R)X'_1 - L'\xi - X'_0)w \right) \partial_\xi \varphi \right) (t, \xi) \, dt \, d\xi \\ \qquad + \int_0^1 L(0) u^{\text{init}}(0, L(0)\xi) \varphi(0, \xi) \, d\xi + \int_0^1 (a - bw(t, 0)) \varphi(\cdot, 0) \, dt = 0; $$

* $X'_1(t) = -\alpha_1 + \beta_1 w(t, 1)$ and $X'_0(t) = \alpha_0 - \beta_0 w(t, 0) + (1-R)X'_1(t)$;

* and $u(t, x) = w(t, \frac{x - X_0(t)}{L(t)})$.

In order to address the convergence of the scheme, let us now explain how we define some approximates solutions based on the discrete data $((u_\ell^n)_{0 \le i \le I+1}, X_0^n, X_1^n)_{0 \le n \le N}$. For the boundaries, we set

$$ (38) \quad X_{\ell,h,\Delta t}(t) = X_\ell^n \quad \text{if } t \in (t^{n-1}, t^n], \quad \text{and} \quad X_{\ell,h,\Delta t}(0) = X_\ell^0, \quad \ell \in \{0, 1\}, $$

so that we can also define $L_{h,\Delta t} = X_{1,h,\Delta t} - X_{0,h,\Delta t}$. Then, it is possible to reconstruct the approximate densities on $\bigcup_{t \in [0,T]} (\{t\} \times [X_{0,h,\Delta t}(t), X_{1,h,\Delta t}(t)])$ by setting

$$ u_{h,\Delta t}(t,x) = \begin{cases} u_i^n & \text{if } (t,x) \in (t^{n-1}, t^n] \times (x_{i-\frac{1}{2}}^n, x_{i+\frac{1}{2}}^n) \\ u_i^0 & \text{if } t=0 \text{ and } x \in (x_{i-\frac{1}{2}}^0, x_{i+\frac{1}{2}}^0) \\ u_\ell^n & \text{if } x = X_\ell^n \text{ and } t \in (t^{n-1}, t^n], \text{ for } \ell \in \{0, 1\}. \end{cases} $$

This reconstruction is useful to show simultaneously the evolution of the oxide layer and the evolution of the density profile inside the layer, as on Figure 2. However, in order to compare different approximate solutions –when studying the long time behavior of the scheme, or when studying its
convergence–, it seems relevant to reconstruct the approximate densities on the underlying fixed domain $[0, T] \times [0, 1]$. This approximate solution will be denoted by $w_{h,\Delta t}$ and defined by

$$ (39) \qquad w_{h,\Delta t}(t, \xi) = \begin{cases} u_i^n & \text{if } (t, \xi) \in (t^{n-1}, t^n] \times (\xi_{i-\frac{1}{2}}, \xi_{i+\frac{1}{2}}) \\ u_i^0 & \text{if } t=0 \text{ and } \xi \in (\xi_{i-\frac{1}{2}}, \xi_{i+\frac{1}{2}}) \\ u_i^n & \text{if } \xi = \ell \text{ and } t \in (t^{n-1}, t^n], \text{ for } \ell \in \{0, 1\}. \end{cases} $$

We will establish the convergence of the scheme by studying the convergence of $(w_{h,\Delta t}, X_{0,h,\Delta t}, X_{1,h,\Delta t})$ when $h$ and $\Delta t$ tend to 0. As a by-product of the convergence of the scheme, we will obtain the existence of a weak solution to (1) in the sense of Definition 4 thanks to Remark 4.1.

In the spirit of [20], we introduce the reconstruction operator $\Pi_M$ corresponding to a mesh $M = (\xi_{i+\frac{1}{2}})_{0 \le i \le I}$, mapping $\mathbb{R}^{I+2}$ to bounded functions defined on $[0, 1]$. We avoid here the explicit introduction of the reconstruction operator in the notation, but rather denote by $z_h = \Pi_M(z)$ the image of $z = (z_i)_{0 \le i \le I+1} \in \mathbb{R}^{I+2}$, defined by

$$ \Pi_M z(\xi) = z_i \text{ for } \xi \in (\xi_{i-\frac{1}{2}}, \xi_{i+\frac{1}{2}}), \quad \Pi_M z(0) = z_0, \text{ and } \Pi_M z(1) = z_{I+1}. $$

We define the finite dimensional vector space $S_M = \text{Range}(\Pi_M)$, equipped with discrete $H^1(0, 1)$ and $(H^1)'(0, 1)$ norms. More precisely, for all $z_h = \Pi_M z \in S_M$, we set

$$ \begin{aligned} \|z_h\|_1 &= \left( \sum_{i=1}^{I+1} \frac{(z_i - z_{i-1})^2}{h_{i-\frac{1}{2}}} + z_0^2 \right)^{\frac{1}{2}}, \\ \|z_h\|_{-1} &= \sup \left\{ \int_0^1 z_h \eta_h / \eta_h \in S_M, \|\eta_h\|_1 \le 1 \right\}, \end{aligned} $$

We recall the discrete Poincaré inequality satisfied in this framework ([12]):

$$ (40) \qquad \|z_h\|_{L^2(0,1)} \le \sqrt{2} \|z_h\|_1, \quad z_h \in S_M. $$

For time iterated discrete vectors $z = (z^n)_{0 \le n \le N} \in (\mathbb{R}^{I+2})^{N+1}$, we define the linear reconstruction operator $\Pi_{M,\Delta t}$ mapping $(\mathbb{R}^{I+2})^{N+1}$ into the bounded functions on $[0, T] \times [0, 1]$ by setting

$$ \Pi_{M,\Delta t} z(t, \cdot) = \Pi_M z^n \quad \text{if } t \in (t^{n-1}, t^n], \quad t \in (0, T], \text{ and } \Pi_M z^0(\xi) = z_h^0. $$

The space $S_{M,\Delta t} = \text{Range}(\Pi_{M,\Delta t})$ is equipped with the discrete $L^2(H^1)$ and $L^2(H^{-1})$ norms

$$ \|\Pi_{M,\Delta t} z\|_{2,1} = \left( \sum_{n=1}^{N} \Delta t \|\Pi_M z^n\|_1^2 \right)^{\frac{1}{2}}, \quad \|\Pi_{M,\Delta t} z\|_{2,-1} = \left( \sum_{n=1}^{N} \Delta t \|\Pi_M z^n\|_{-1}^2 \right)^{\frac{1}{2}}. $$

Denoting by $u = (u^n)_{0 \le n \le N} \in (\mathbb{R}^{I+2})^{N+1}$ with $u^n = (u_i^n)_{0 \le i \le I+1}$, we have that $w_{h,\Delta t}$ defined by (39) verifies $w_{h,\Delta t} = \Pi_{M,\Delta t}(u)$ with $w_h^n = \Pi_M(u^n)$.

We also define the approximate time and space derivatives of $w_{h,\Delta t}$, respectively $\partial_{t,\Delta t}w_{h,\Delta t}$ and $\partial_{\xi,h}w_{h,\Delta t}$, by

$$ (41) \qquad \partial_{t,\Delta t} w_{h,\Delta t} = \begin{cases} \frac{w_{h,\Delta t}(t)-w_{h,\Delta t}(t-\Delta t)}{\Delta t} & \text{for } t \in (\Delta t, T], \\ \frac{w_h^1-w_h^0}{\Delta t} & \text{for } t \in [0, \Delta t]. \end{cases} $$

and

$$ (42) \qquad \partial_{\xi,h} w_{h,\Delta t}(t, \xi) = \frac{w_{i+1}^n - w_i^n}{h_{i+\frac{1}{2}}} \quad \text{if } (t, \xi) \in (t^{n-1}, t^n] \times (\xi_i, \xi_{i+1}). $$
Note that $\partial_{t,\Delta t} w_{h,\Delta t}$ belongs to $S_{M,\Delta t}$, but that $\partial_{\xi,h} w_{h,\Delta t}$ does not.

Finally, we define the approximate flux $\mathcal{F}_{h,\Delta t} \in L^\infty((0,T) \times (0,1))$ by

$$ (43) \qquad \mathcal{F}_{h,\Delta t}(t, \xi) = \mathcal{F}_{i+\frac{1}{2}}^n \quad \text{if } (t, \xi) \in (t^{n-1}, t^n] \times (\xi_i, \xi_{i+1}). $$

4.2. **Additional estimates.** In this section, we derive some additional energy estimates on the approximate solution on a time interval $(0, T)$ with the time horizon $T$ being such that there exists some uniform bounds $L^\# \ge L_\# > 0$ for the domain size, i.e. $L^\# \ge L^n \ge L_\#$ for all $n \in \{0, \dots, N\}$. The existence of such a positive $T$ is granted by Proposition 3.6.

**Proposition 4.2.** Assume that Assumption (12) holds. Let $T > 0$ be a finite time horizon such that $L^\# > L^n \ge L_\# > 0$ for all $n \in \{0, \dots, N\}$, then there exists $C_1 > 0$ depending only on $L^\#$, $m$, $M$, on the parameters $a$, $b$, $\alpha_0$, $\beta_0$, $\alpha_1$, $\beta_1$, $R$ and $T$ (but neither on the mesh $M$ nor on $\Delta t$) such that $\|w_{h,\Delta t}\|_{2,1}^2 \le C_1$.

*Proof.* Proposition 3.2 gives

$$ (44) \qquad \frac{\mathcal{H}_{\phi}^{\text{tot},n} - \mathcal{H}_{\phi}^{\text{tot},n-1}}{\Delta t} + \mathcal{D}_{\phi}^{\text{bulk},n} \le -\mathcal{D}_{\phi}^{\text{bound},n} \le 0, \quad 1 \le n \le N, $$

with $\mathcal{D}_{\phi}^{\text{bulk},n}$ and $\mathcal{D}_{\phi}^{\text{bound},n}$ being defined by (25) and (27) respectively. For the specific choice $\phi(r) = \frac{1}{2}r^2$, we have

$$ \mathcal{D}_{\phi}^{\text{bulk},n} = \sum_{i=0}^{I} \frac{(u_{i+1}^n - u_i^n)^2}{L^n h_{i+\frac{1}{2}}} \left( B \left( -L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n \right) \theta_{i+\frac{1}{2}}^n + B \left( L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n \right) (1 - \theta_{i+\frac{1}{2}}^n) \right). $$

Since $B$ is decreasing, one has $B(r) \ge B(|r|)$. Then since $\theta_{i+\frac{1}{2}}^n \in [0, 1]$ and since $L^n \le L^\#$, we obtain that

$$ \mathcal{D}_{\phi}^{\text{bulk},n} \ge \frac{1}{L^{\#}} \sum_{i=0}^{I} \frac{(u_{i+1}^n - u_i^n)^2}{h_{i+\frac{1}{2}}} B \left( L^n h_{i+\frac{1}{2}} |v_{i+\frac{1}{2}}^n| \right). $$

Thanks to the bounds on $v_{i+\frac{1}{2}}^n$ stated in Corollary 3.4 combined with the upper bound $L^\#$ on $L^n$ and with $h_{i+\frac{1}{2}} \le 1$, we get some uniform upper bound on $L^n h_{i+\frac{1}{2}} |v_{i+\frac{1}{2}}^n|$, whence a uniform lower bound $\delta > 0$ for $B \left( L^n h_{i+\frac{1}{2}} |v_{i+\frac{1}{2}}^n| \right)$. As a consequence, estimate (44) yields

$$ \frac{\delta}{L^{\#}} \|w_h^n\|_1^2 \le \mathcal{D}_{\phi}^{\text{bulk},n} + \frac{\delta}{L^{\#}} (u_0^n)^2 \le -\frac{\mathcal{H}_{\phi}^{\text{tot},n} - \mathcal{H}_{\phi}^{\text{tot},n-1}}{\Delta t} + \frac{\delta}{L^{\#}} M^2. $$

Then, multiplying by $\Delta t$ and summing over $n \in \{1, \dots, N\}$ gives

$$ (45) \qquad \|w_{h,\Delta t}\|_{2,1}^2 \le L^{\#} \frac{\mathcal{H}_{\phi}^{\text{tot},0} - \mathcal{H}_{\phi}^{\text{tot},N}}{\delta} + M^2 T. $$

We infer from the definition (19) of the discrete total free energy that

$$ (46) \qquad \mathcal{H}_{\phi}^{\text{tot},0} = \mathcal{H}_{\phi}^0 \le \mathcal{H}_{\phi}(0) \le \frac{L^{\#} M^2}{2}. $$

On the other hand, since $\mathcal{H}_{\phi}^N \ge 0$, one gets from (19) combined with the uniform upper bound (29) on $u_i^n$ that

$$ (47) \qquad \mathcal{H}_{\phi}^{\text{tot},N} \ge -T \left( \frac{1}{2} \left( \frac{\alpha_0}{\beta_0} \right)^2 \left( M - \frac{\alpha_0}{\beta_0} \right) + a \left( M - \frac{a}{b} \right) + \frac{R}{2} \left( \frac{\alpha_1}{\beta_1} \right)^2 \left( M - \frac{\alpha_1}{\beta_1} \right) \right). $$
Taking (46) and (47) into account in (45) concludes the proof.

□

Our next statement is about the control in the discrete $L^2((H^1)')$ norm of the discrete time derivative $\partial_{t,\Delta t}w_{h,\Delta t}$ introduced in (41).

**Proposition 4.3.** Under the assumptions of Proposition 4.2, there exist $C_2 > 0$ and $C_3 > 0$ depending only on $L_b, L^\sharp, m, M$, on the parameters $a, b, \alpha_0, \beta_0, \alpha_1, \beta_1, R$ and $T$, but not on the discretization, such that

$$ (48) \qquad \|\mathcal{F}_{h,\Delta t}\|_{L^2((0,T)\times(0,1))} \le C_2 $$

$$ (49) \qquad \|\partial_{t,\Delta t} w_{h,\Delta t}\|_{2,-1} \le C_3. $$

*Proof.* Let $\eta_h \in S_M$ with $\eta_h(\xi) = \eta_i$ for $\xi \in (\xi_{i-\frac{1}{2}}, \xi_{i+\frac{1}{2}})$ be such that $\|\eta_h\|_1 \le 1$, then for $t \in (t^{n-1}, t^n]$, $1 \le n \le N$, one has

$$ E^n = \int_0^1 \partial_{t,\Delta t} w_{h,\Delta t}(t, \xi) \eta_h(\xi) d\xi = \sum_{i=1}^I h_i \frac{u_i^n - u_i^{n-1}}{\Delta t} \eta_i. $$

Rewriting

$$ \frac{u_i^n - u_i^{n-1}}{\Delta t} = \frac{1}{L^n} \left( \frac{L^n u_i^n - L^{n-1} u_i^{n-1}}{\Delta t} - u_i^{n-1} \frac{L^n - L^{n-1}}{\Delta t} \right), $$

and using (6a), one gets that

$$ (50) \qquad E^n = -\frac{1}{L^n} \sum_{i=1}^{I} \eta_i (\mathcal{F}_{i+\frac{1}{2}}^n - \mathcal{F}_{i-\frac{1}{2}}^n) - \frac{1}{L^n} \frac{L^n - L^{n-1}}{\Delta t} \sum_{i=1}^{I} \eta_i u_i^{n-1} h_i. $$

As a direct consequence of (31), there exists some $C_4 > 0$ depending only on the data of the continuous problem such that $|\delta^n[L]| \le C_4$. The Poincaré inequality (40), the lower bound on $L^n$ and the upper bound on $u_i^{n-1}$ then yield

$$ \left| \frac{1}{L^n} \frac{L^n - L^{n-1}}{\Delta t} \sum_{i=1}^{I} \eta_i u_i^{n-1} h_i \right| \le \frac{C_4}{L_b} \|w_h^n\|_{L^2(0,1)} \|\eta_h\|_{L^2(0,1)} \le \frac{C_4}{L_b} \sqrt{2M}. $$

Besides, applying a summation by parts to the first term in the right-hand-side of (50) provides

$$ -\frac{1}{L^n} \sum_{i=1}^{I} \eta_i (\mathcal{F}_{i+\frac{1}{2}}^n - \mathcal{F}_{i-\frac{1}{2}}^n) = \frac{1}{L^n} \sum_{i=0}^{I} \mathcal{F}_{i+\frac{1}{2}}^n (\eta_{i+1} - \eta_i) + \frac{1}{L^n} \mathcal{F}_{\frac{1}{2}}^n \eta_0. $$

Since $|\mathcal{F}_{\frac{1}{2}}^n| \le \min(a - bm, bM - a)$, and since $|\eta_0| \le \|\eta_h\|_1 \le 1$, one obtains that

$$ \left| \frac{1}{L^n} \mathcal{F}_{\frac{1}{2}}^n \eta_0 \right| \le \frac{\min(a - bm, bM - a)}{L_b}. $$

On the other hand, applying the Cauchy-Schwarz inequality, one gets that

$$ \frac{1}{L^n} \sum_{i=0}^{I} \mathcal{F}_{i+\frac{1}{2}}^n (\eta_{i+1} - \eta_i) \le \frac{1}{L_b} \left( \sum_{i=0}^{I} h_{i+\frac{1}{2}} |\mathcal{F}_{i+\frac{1}{2}}^n|^2 \right)^{\frac{1}{2}}. $$

Thanks to relation (9), the Scharfetter-Gummel fluxes rewrite

$$ (51) \qquad \mathcal{F}_{i+\frac{1}{2}}^n = v_{i+\frac{1}{2}}^n \frac{u_i^n + u_{i+1}^n}{2} + \frac{1}{2L^n} \left( B(-L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n) + B(L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n) \right) \frac{u_i^n - u_{i+1}^n}{h_{i+\frac{1}{2}}}, $$
The velocities $v_{i+\frac{1}{2}}^n$ are uniformly bounded owing to (32). As highlighted in the proof of Proposition 4.2, this yields a uniform bound on the quantity

$$ \frac{1}{2L^n} \left( B(-L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n) + B(L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n) \right) $$

as well. Using the elementary inequality $(p+q)^2 \le 2(p^2+q^2)$, we therefore obtain that

$$ \sum_{i=0}^{I} h_{i+\frac{1}{2}} \left| \mathcal{F}_{i+\frac{1}{2}}^{n} \right|^{2} \le C_{5} + C_{6} \| w_{h}^{n} \|_{1}^{2} $$

for some quantities $C_5$ and $C_6$ not depending on the discretization parameters, but only on the data prescribed in the statement of Proposition 4.3. Summing the above estimate over $n \in \{1, \dots, N\}$ and using Proposition 4.2 gives exactly (48).

Summing up the above calculations, we end up with the existence of quantities $C_7$ and $C_8$ depending only on the prescribed data such that $|E^n|^2 \le C_7 + C_8 \|w_h^n\|_1^2$ whatever the choice of $\eta_h \in S_M$ provided that $\|\eta_h\|_1 \le 1$. Therefore, this implies that

$$ \left\| \frac{w_h^n - w_h^{n-1}}{\Delta t} \right\|_{-1}^2 \le C_7 + C_8 \|w_h^n\|_1^2. $$

Multiplying by $\Delta t$, summing over $n$, and using Proposition 4.2 proves (49). $\square$

### 4.3. Compactness results.
In this part is presented the main compactness result for the solution to the numerical scheme (6). Before focusing on the density profile, let us start by the interfaces $X_{h,\Delta t} = (X_{0,h,\Delta t}, X_{1,h,\Delta t})$ and their approximate derivatives

$$ \partial_{t,\Delta t} X_{\ell,h,\Delta t}(t) = \begin{cases} \frac{X_{\ell,h,\Delta t}(t) - X_{\ell,h,\Delta t}(t-\Delta t)}{\Delta t} & \text{for } t \in (\Delta t, T], \\ \frac{X_{\ell}^1 - X_{\ell}^0}{\Delta t} & \text{for } t \in [0, \Delta t], \end{cases} \quad \ell \in \{0, 1\}. $$

**Lemma 4.4.** Under the assumptions of Proposition 4.2, there exist some Lipschitz continuous functions $X_0$ and $X_1$ defined on $[0, T]$, with $X_0(0) = 0$, $X_1(0) = L^0$ and with $L(t) = X_1(t) - X_0(t) \in [L_\nu, L^\ddag]$ for all $t \in [0, T]$ such that, up to a subsequence

$$ (52) \qquad X_{\ell,h,\Delta t} \xrightarrow{h, \Delta t \to 0} X_{\ell} \quad \text{uniformly on } [0, T], \quad \ell \in \{0, 1\}, $$

while

$$ (53) \qquad \partial_{t,\Delta t} X_{\ell,h,\Delta t} \xrightarrow{h,\Delta t \to 0} X'_{\ell} \quad \text{in the } L^{\infty}(0,T)\text{-weak-\star sense.} $$

**Proof.** The uniform bounds on $\Delta X_\ell^n$ yield uniform bounds on $\partial_{t,\Delta t} X_{\ell,h,\Delta t}$, thus there exist $g_\ell$ in $L^\infty(0,T)$, $\ell \in \{0,1\}$, such that

$$ \partial_{t,\Delta t} X_{\ell,h,\Delta t} \xrightarrow{h,\Delta t \to 0} g_\ell \quad \text{in the } L^\infty(0,T)\text{-weak-\star sense.} $$

For $\tau > 0$ and $t \in [0, T - \tau)$, one furthermore has

$$ |X_{\ell,h,\Delta t}(t+\tau) - X_{\ell,h,\Delta t}| \le C(\tau + \Delta t) $$

with $C$ being a uniform upper bound for $\Delta X_\ell^n$, $n \ge 0$. We can then apply the refined version of the Ascoli-Arzelà theorem [2, Proposition 3.3.1] to infer the desired uniform convergence. The preservation of the bounds on $L$ follows. Showing that $g_\ell = X'_\ell$ in the sense of distributions in then classical. $\square$
Let us now consider the compactness of the approximate density profiles on the fixed domain $w_{h,\Delta t}$.

**Lemma 4.5.** Under the assumptions of Proposition 4.2, there exists a function $w \in L^\infty((0,T) \times (0,1)) \cap L^2(0,T); H^1(0,1))$ such that, up to a subsequence, $w_{h,\Delta t} \to w$ (strongly) in $L^2((0,T) \times (0,1))$ as $h, \Delta t \to 0$ and $\partial_{\xi,h}w_{h,\Delta t} \to \partial_{\xi}w$ weakly in $L^2((0,T) \times (0,1))$. Moreover, we have the strong convergence of the traces

$$ (54) \qquad w_{h,\Delta t}(t, 0) \underset{h,\Delta t \to 0}{\longrightarrow} w(t, 0) \quad \text{and} \quad w_{h,\Delta t}(t, 1) \underset{h,\Delta t \to 0}{\longrightarrow} w(t, 1) \quad \text{a.e. in } (0,T), $$

and initial profile

$$ (55) \qquad w_{h,\Delta t}(0, \xi) \underset{h \to 0}{\longrightarrow} w^{\text{init}}(\xi) = u^{\text{init}}(\xi L^0) \quad \text{for a.e. } \xi \in (0,1). $$

*Proof.* As $m \le w_{h,\Delta t} \le M$ a.e. in $(0,T) \times (0,1)$, there exists $w$ measurable with $m \le w \le M$ such that $w_{h,\Delta t}$ converges, up to a subsequence, towards $w$ in the $L^\infty(0,T) \times (0,1)$ weak-⋆ sense. With Propositions 4.2 and 4.3 at hand, one can directly apply the discrete Aubin-Simon lemma [23, Theorem 3.4], and get that, up to the extraction of another subsequence, $w_{h,\Delta t} \to w$ in $L^2((0,T); L^2(0,1))$ as $h, \Delta t \to 0$. Moreover, we deduce from Proposition 4.2 that $\partial_{\xi,h}w_{h,\Delta t}$ is bounded in $L^2((0,T) \times (0,1))$, hence its weak in $L^2$ convergence towards some limit $f$. A straightforward adaptation in [21, Lemma 8.2] shows that the limiting profile $w$ belongs to $L^2((0,T); H^1(0,1))$. The identification $f = \partial_{\xi}w$ is classical.

Finally, the convergence of the traces (54) can be proven by adapting the strategy of [8, Section 4.2] to our simpler one-dimensional setting, while the convergence property (55) is nothing but the classical consistence result for the approximation (5). □

**4.4. Identification of the limit as a weak solution.** The goal of this section is to establish our last theoretical result, which by the way grants the existence of a weak solution up to a minimal time thanks to Proposition 3.6. As a preliminary, let us show some compactness and consistency result on the approximate fluxes following ideas from [13, Lemma 4.4].

**Lemma 4.6.** Under the assumptions of Proposition 4.2, $\mathcal{F}_{h,\Delta t}$ defined by (43) converges, up to a subsequence, weakly in $L^2((0,T) \times (0,1))$ towards

$$ (56) \qquad \mathcal{F} = vw - \frac{1}{L} \partial_{\xi}w \quad \text{with} \quad v(t, \xi) = (1-R)X_1'(t) - L'(t)\xi - X_0'(t). $$

*Proof.* As a consequence of Lemma 4.4, the discrete velocity field

$$ v_{h,\Delta t}(t, \xi) = v_{i+\frac{1}{2}}^n \quad \text{if } (t, \xi) \in (t^{n-1}, t^n] \times (\xi_i, \xi_{i+1}) $$

converges in the $L^\infty((0,T) \times (0,1))$ weak-⋆ sense towards $v$ defined in (56).

Starting from expression (51) for the flux $\mathcal{F}_{i+\frac{1}{2}}^n$ and using $B(0) = 1$, one gets that

$$ (57) \qquad \mathcal{F}_{h,\Delta t} = v_{h,\Delta t} \bar{w}_{h,\Delta t} - \frac{1}{L_{h,\Delta t}} \partial_{\xi,h} w_{h,\Delta t} + r_{h,\Delta t}, $$

where we have set

$$ r_{h,\Delta t}(t, \xi) = \frac{u_i^n - u_{i+1}^n}{2L^n h_{i+\frac{1}{2}}} \left( B(-L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n) + B(L^n h_{i+\frac{1}{2}} v_{i+\frac{1}{2}}^n) - 2B(0) \right) $$

and $\bar{w}_{h,\Delta t}(t, \xi) = \frac{u_i^n + u_{i+1}^n}{2}$ for $(t, \xi) \in (t^{n-1}, t^n] \times (\xi_i, \xi_{i+1})$. Since $B$ is 1-Lipschitz, we obtain that

$$ |r_{h,\Delta t}| \le h |v_{h,\Delta t}| |\partial_{\xi,h} w_{h,\Delta t}| \underset{h,\Delta t \to 0}{\longrightarrow} 0 \quad \text{in } L^2((0,T) \times (0,1)) $$
thanks to the uniform boundedness of $v_{h,\Delta t}$ in $L^\infty(0,T) \times (0,1))$ and of $\partial_{\xi,h}w_{h,\Delta t}$ in $L^2((0,T) \times (0,1))$, cf. Proposition 4.2. Moreover, one readily checks that

$$|\bar{w}_{h,\Delta t}(t, \xi) - w_{h,\Delta t}(t, \xi)| \leq h|\partial_{\xi,h}w_{h,\Delta t}|,$$

so that, thanks to Lemma 4.5, $\bar{w}_{h,\Delta t}(t, \xi)$ also converges strongly in $L^2$ towards $w$. So, using Lemmas 4.4 ad 4.5, we can pass to the limit $h, \Delta t \to 0$ in (57) to recover (56). $\square$

We are now in position to prove the main result of this section.

**Theorem 4.7.** Let $X = (X_0, X_1)$ be as in Lemma 4.4, let $w$ be defined as in Lemma 4.5 and $u$ be defined by $u(t,x) = w(t, \frac{x-X_0(t)}{X_1(t)-X_0(t)})$, then $(u, X_0, X_1)$ is a weak solution in the sense of Definition 4.

*Proof.* Relations (6d)–(6e) for the motion of the free boundaries rewrite

$$\begin{aligned} \partial_{t,\Delta t} X_{1,h,\Delta t} &= -\alpha_1 + \beta_1 w_{h,\Delta t}(\cdot, 1), \\ \partial_{t,\Delta t} X_{0,h,\Delta t} &= (1-R)\partial_{t,\Delta t} X_{1,h,\Delta t} + \alpha_0 - \beta_0 w_{h,\Delta t}(\cdot, 0). \end{aligned}$$

We can pass to the limit $h, \Delta t \to 0$ in the above relations thanks to Lemmas 4.4 and 4.5, leading to $X_1' = -\alpha_1 + \beta_1 w(\cdot, 1)$ and $X_0' = (1-R)X_1' + \alpha_0 - \beta_0 w(\cdot, 0)$.

Now let $\varphi \in C_c^\infty([0,T) \times [0,1])$. For $i \in \{0, \dots, I+1\}$ and $n \in \{0, \dots, N\}$, we denote by $\varphi_i^n = \varphi(t^n, \xi_i)$, by $\varphi_{h,\Delta t} = \Pi_{M,\Delta t}((\varphi_i^n)_{i,n})$, and by $\partial_{t,\Delta t}\varphi_{h,\Delta t}$ and $\partial_{\xi,h}\varphi_{h,\Delta t}$ the approximate time and space derivatives of $\varphi_{h,\Delta t}$ respectively defined with formulas (41) and (42) with $\varphi$ instead of $w$. Note that $\varphi_i^n = 0$. Because of the regularity of $\varphi$, one readily shows that

$$(58) \qquad \varphi_{h,\Delta t} \underset{h,\Delta t \to 0}{\longrightarrow} \varphi \quad \text{uniformly on } [0,T] \times [0,1],$$

$$(59) \qquad \partial_{t,\Delta t} \varphi_{h,\Delta t} \underset{h,\Delta t \to 0}{\longrightarrow} \partial_t \varphi \quad \text{in the } L^\infty((0,T) \times (0,1))\text{-weak-}^* \text{ sense,}$$

$$(60) \qquad \partial_{\xi,h} \varphi_{h,\Delta t} \underset{h,\Delta t \to 0}{\longrightarrow} \partial_{\xi} \varphi \quad \text{uniformly on } [0,T] \times [0,1].$$

Multiplying the scheme (6a) by $\varphi_i^{n-1}$ and summing over $1 \leq i \leq I$ and $1 \leq n \leq N$ yields

$$(61) \qquad A_{h,\Delta t} + B_{h,\Delta t} = 0,$$

where we have set

$$\begin{aligned} A_{h,\Delta t} &= \sum_{n=1}^{N} \sum_{i=1}^{n} (L^n u_i^n - L^{n-1} u_i^{n-1}) \varphi_i^{n-1} h_i, \\ B_{h,\Delta t} &= \sum_{n=1}^{N} \Delta t \sum_{i=1}^{n} \left( \mathcal{F}_{i+\frac{1}{2}}^{n} - \mathcal{F}_{i-\frac{1}{2}}^{n} \right) \varphi_i^{n-1}. \end{aligned}$$

Applying a summation by parts to $A_{h,\Delta t}$ gives

$$A_{h,\Delta t} = A'_{h,\Delta t} + A''_{h,\Delta t} = \sum_{n=1}^{N} \sum_{i=1}^{n} L^n u_i^n (\varphi_i^{n-1} - \varphi_i^n) h_i - \sum_{i=1}^{n} L^0 u_i^0 \varphi_i^0.$$

Because of the uniform convergence of $L_{h,\Delta t}$ towards $L$, cf. Lemma 4.4, of the strong in $L^2$ convergence of $w_{h,\Delta t}$ towards $w$, cf. Lemma 4.5, and thanks to (59), one can pass to the limit in
the first term in the above right-hand side to obtain

$$A'_{h,\Delta t} = \int_0^T \int_0^1 L_{h,\Delta t} w_{h,\Delta t} \partial_{t,\Delta t} \varphi_{h,\Delta t} \xrightarrow{h,\Delta t \to 0} \int_0^T \int_0^1 L w \partial_t \varphi.$$

For the second contribution, one deduces from (55) and (58) that

$$A''_{h,\Delta t} = \int_0^1 L^0 w_h^0 \varphi_{h,\Delta t}(0, \cdot) \xrightarrow{h \to 0} \int_0^1 L^0 u^{\text{init}}(L^0 \xi) \varphi(0, \xi) d\xi.$$

On the other hand, yet another summation by parts provides

$$B_{h,\Delta t} = B'_{h,\Delta t} + B''_{h,\Delta t} = \sum_{n=1}^{N} \Delta t \sum_{i=1}^{n} \mathcal{F}_{i+\frac{1}{2}}^{n} (\varphi_{i}^{n-1} - \varphi_{i+1}^{n}) - \sum_{n=1}^{N} \mathcal{F}_{\frac{1}{2}}^{n} \varphi_{0}^{n}.$$

The first term in the above right-hand side rewrites

$$B'_{h,\Delta t} = - \int_0^T \int_0^1 \mathcal{F}_{h,\Delta t} \partial_{\xi,h} \varphi_{h,\Delta t} \xrightarrow{h,\Delta t \to 0} \int_0^T \int_0^1 \left( \frac{1}{L} \partial_{\xi} w - vw \right) \partial_{\xi} \varphi$$

owing to Lemma 4.6 and to (60). For the second term, one makes use of the definition (6b) of the flux $\mathcal{F}_{\frac{1}{2}}^n$ to write

$$B''_{h,\Delta t} = \int_0^T (b w_{h,\Delta t}(t, 0) - a) \varphi_{h,\Delta t}(t, 0) dt \xrightarrow{h,\Delta t \to 0} \int_0^T (b w(t, 0) - a) \varphi(t, 0) dt$$

thanks to (58), to (54). We then deduce from (61) and from the above calculations that $(w, X_0, X_1)$ satisfies (37) and all the points of Remark 4.1 are satisfied. It concludes the proof of Theorem 4.7. $\square$

## 5. NUMERICAL RESULTS

In this final section, we present some numerical experiments. The scheme being implicit, its solutions are obtained thanks to Newton's method, with a $10^{-10}$ tolerance between two consecutive iterations, on a uniform mesh. The different test cases, corresponding to different parameters, are presented in Table 1. The numerical results are obtained with $I = 100$ cells, and the time step is $\Delta t = 10^{-2}$, until different final times depending on the test case.

Test case 1 satisfies the condition (12), and one observes on Figure 2 the convergence of the approximate solution $u_{h,\Delta t}$ towards the travelling wave profile $\hat{u}$. Condition (12) is satisfied neither for test case 2, for which the oxide layer disappears in finite time, nor for test case 3 where the oxide layer grows indefinitely.

<table><caption>TABLE 1. Parameters of the test cases.</caption><thead><tr><th></th><th>a</th><th>b</th><th>α<sub>0</sub></th><th>β<sub>0</sub></th><th>α<sub>1</sub></th><th>β<sub>1</sub></th><th>R</th><th>u<sup>init</sup></th><th>L<sup>0</sup></th><th>T</th><th>(12)</th></tr></thead><tbody><tr><td>Test case 1</td><td>1</td><td>1</td><td>1.5</td><td>1</td><td>0.5</td><td>4</td><td>2</td><td><sup>a</sup>/<sub>b</sub>e<sup>-Rcx</sup> + 2</td><td>1</td><td>20</td><td>✓</td></tr><tr><td>Test case 2</td><td>1.75</td><td>1</td><td>5</td><td>2</td><td>5</td><td>2</td><td>2</td><td><sup>a</sup>/<sub>b</sub>e<sup>-Rcx</sup> + 2</td><td>1</td><td>3.5</td><td>✗</td></tr><tr><td>Test case 3</td><td>1</td><td>1</td><td>4</td><td>1</td><td>3</td><td>1.5</td><td>2</td><td><sup>a</sup>/<sub>b</sub>e<sup>-Rcx</sup> + 2</td><td>1</td><td>10</td><td>✗</td></tr></tbody></table>

In order to compare the travelling wave profile and the other profile from test case 1, both are rescaled in $[0, 1]$, so it is possible compute a discrete distance between the approximate solution $u$ and the travelling wave profile $\hat{u}$. Precisely, after the rescaling, the distance between the two profiles at each time step $n$ is measured through the quantity $d^n = L^n \sum_{i=1}^I h_i (u_i^n - \hat{u}(\xi_i))^2$. The
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Profiles of the approximate solution for test cases 1 (left), 2 (center) and 3 (right).</figcaption></figure>

evolution of this distance along time is plotted at the left of Figure 3 in a semilogarithmic scale, to exhibit the exponential convergence of the (approximate) solution towards a travelling wave as time goes to $+\infty$. The evolution of the distance between the approximate width of the oxide $L^n$ and the one of the travelling wave $\hat{L}$ is also plotted, at the right of Figure 3. This behavior is in good accordance with the results obtained in [4] and [5] for more complex corrosion models. Observe on Figure 3 that the computation ultimately reaches machine error. This is a consequence of the exactness of our scheme for the travelling wave.

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Convergence towards a travelling wave profile of the concentration $u$ (left) and the width of the oxide $L$ (right).</figcaption></figure>

Finally, we compute approximate solutions $(w_{h_k, \Delta t_k}, X_{0,h_k, \Delta t_k}, X_{1,h_k, \Delta t_k})$ on different meshes of $I_k = 50 \times 2^k$ cells and $N_k = 10 \times 4^k$ time steps for $k \in \{0, \dots, 5\}$ for test-case 1. As we already know that the scheme is exact for the travelling wave, which has been shown on Figure 3 to be the long-time limit the scheme, we consider a rather small final time $T = 0.2$ to focus on the transient regime. The reference solution $(w_{\text{ref}}, X_{0,\text{ref}}, X_{1,\text{ref}})$ is the one on the finest mesh, i.e. for $k = 5$. Introducing the orthogonal projection $P_{M_k, \Delta t_k} : L^2((0,T) \times (0,1)) \to S_{M_k, \Delta t_k}$, then the error and the associated convergence rate for the profile $w_{h_k, \Delta t_k}$ is estimated by

$$ \begin{aligned} \text{err}_{w,k} &= \|w_{h_k, \Delta t_k} - P_{M_k, \Delta t_k} w_{\text{ref}}\|_{L^2((0,T) \times (0,1))}, \\ \text{rate}_{w,k} &= \frac{\log(\text{err}_{w,k}) - \log(\text{err}_{w,k-1})}{\log h_k - \log h_{k-1}}. \end{aligned} $$

Similarly, denoting by $P_{\Delta t_k}$ the orthogonal projection of $L^2(0,T)$ on the set of piecewise constant functions on time intervals of constant length $\Delta t_k$, we denote by

$$ \text{err}_{\ell,k} = \|X_{\ell,h_k, \Delta t_k} - P_{\Delta t_k} X_{\ell, \text{ref}}\|_{L^{\infty}(0,T)}, \quad \ell \in \{0,1\}. $$
As the interface location are functions of time only, we choose to compute the order w.r.t. the time discretization, i.e.,

$$ \text{rate}_{\ell,k} = \frac{\log(\text{err}_{\ell,k}) - \log(\text{err}_{\ell,k-1})}{\log \Delta t_k - \log \Delta t_{k-1}}, \quad \ell \in \{0, 1\}, $$

which is in our case twice smaller than the one w.r.t. the space discretization. One observes in Table 2 that the scheme is second order accurate in space and first order accurate in time, as expected from the choice of Scharfetter-Gummel fluxes and of the backward Euler time integration scheme.

<table><caption>TABLE 2. Convergence of the scheme under grid refinement</caption><thead><tr><th>k</th><th>h<sub>k</sub></th><th>Δt<sub>k</sub></th><th>err<sub>w,k</sub></th><th>rate<sub>w,k</sub></th><th>err<sub>0,k</sub></th><th>rate<sub>0,k</sub></th><th>err<sub>1,k</sub></th><th>rate<sub>1,k</sub></th></tr></thead><tbody><tr><td>0</td><td>2.0e-2</td><td>2.0e-2</td><td>3.36e-3</td><td>-</td><td>7.85e-3</td><td>-</td><td>5.40e-3</td><td>-</td></tr><tr><td>1</td><td>1.0e-2</td><td>5.0e-3</td><td>8.90e-4</td><td>1.91</td><td>2.22e-3</td><td>0.90</td><td>1.86e-3</td><td>0.76</td></tr><tr><td>2</td><td>5.0e-3</td><td>1.25e-3</td><td>1.96e-4</td><td>2.17</td><td>6.17e-4</td><td>0.92</td><td>5.4e-4</td><td>0.89</td></tr><tr><td>3</td><td>2.5e-3</td><td>3.125e-4</td><td>3.77e-5</td><td>2.38</td><td>1.56e-4</td><td>0.99</td><td>1.4e-4</td><td>0.97</td></tr><tr><td>4</td><td>1.25e-3</td><td>7.81e-5</td><td>5.71e-6</td><td>2.72</td><td>3.18e-5</td><td>1.14</td><td>2.89e-5</td><td>1.14</td></tr></tbody></table>

## 6. CONCLUSION AND PROSPECTS

The work conducted in this paper on a simplified model for corrosion, with moving boundaries, presents both theoretical and numerical results providing rigorous foundations in this simplified setting to phenomena observed numerically on more complex models, especially concerning the crucial role played by travelling waves when such profiles exist.

The numerical strategy employed in the code CALIPSO [5], which also enters the framework of ALE finite volumes, and which also makes use of Scharfetter-Gummel fluxes, is by the way validated. Indeed, the existence of a travelling wave profile as presented in Proposition 2.1, under some explicit assumption on the kinetic parameters, joins the work that is done in [7] and [14]. Moreover, in this one species case, the expression of the travelling wave is explicit, and this motivates the choice of Scharfetter-Gummel fluxes for the discretization employed in [5]. This kind of fluxes guarantees the preservation of the travelling wave profile by the numerical scheme, and is compatible with the rich structure of the continuous problem, which admits an infinite number of free energies.

Beyond the expected convergence orders -2 in space and 1 in time-, the numerical results show the convergence in the long-time regime towards the unique travelling wave. The theoretical analysis of this asymptotic behavior will be conducted in future works. We will also extend our results to richer models for corrosion inspired from [4], including multiple species and self-consistent electric drift.

## ACKNOWLEDGEMENTS

The authors acknowledge partial support from the European Union's Horizon 2020 research and innovation programme under grant agreement No 847593 (WP DONUT) and from Labex CEMPI (ANR-11-LABX-0007-01).
REFERENCES

[1] C. Alamichel and N. Meunier. Existence of traveling wave for a coupled incompressible Darcy's free boundary model with undercooling effect and surface tension. arXiv:2501.04576, 2025.

[2] L. Ambrosio, N. Gigli, and G. Savaré. *Gradient flows in metric spaces and in the space of probability measures.* Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, second edition, 2008.

[3] A. Arnold, P. Markowich, G. Toscani, and A. Unterreiter. On convex Sobolev inequalities and the rate of convergence to equilibrium for Fokker-Planck type equations. *Comm. Partial Differential Equations*, 26(1-2):43–100, 2001.

[4] C. Bataillon, F. Bouchon, C. Chainais-Hillairet, C. Desgranges, E. Hoarau, F. Martin, S. Perrin, M. Turpin, and J. Talandier. Corrosion modelling of iron based alloy in nuclear waste repository. *Electrochim. Acta*, 55(15):4451–4467, 2010.

[5] C. Bataillon, F. Bouchon, C. Chainais-Hillairet, J. Fuhrmann, E. Hoarau, and R. Touzani. Numerical methods for the simulation of a corrosion model with moving oxide layer. *J. Comput. Phys.*, 231(18), 2012.

[6] F. Bolley and J. A. Carrillo. Nonlinear diffusion: geodesic convexity is equivalent to Wasserstein contraction. *Comm. Partial Differential Equations*, 39(10):1860–1869, 2014.

[7] M. Breden, C. Chainais-Hillairet, and A. Zurek. Existence of traveling wave solutions for the Diffusion Poisson Coupled Model: a computer-assisted proof. *ESAIM Math. Model. Numer. Anal.*, 55(4):1669–1697, 2021.

[8] K. Brenner, C. Cancès, and D. Hilhorst. Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure. *Comput. Geosci.*, 17:573–597, 2013.

[9] C. Cancès, J. Cauvin-Vila, C. Chainais-Hillairet, and V. Ehrlacher. Cross-diffusion systems coupled via a moving interface. *Interfaces Free Bound.*, 2024. Online first.

[10] C. Cancès, C. Chainais-Hillairet, B. Merlet, F. Raimondi, and J. Venel. Mathematical analysis of a thermodynamically consistent reduced model for iron corrosion. *Z. Angew. Math. Phys.*, 74(96), 2023.

[11] P. Cermelli and M. Jabbour. Multispecies epitaxial growth on vicinal surfaces with chemical reactions and diffusion. *Proc. R. Soc. A*, 461(2063):3483–3504, 2005.

[12] C. Chainais-Hillairet and C. Bataillon. Mathematical and numerical study of a corrosion model. *Numer. Math.*, 110:1–25, 2008.

[13] C. Chainais-Hillairet and J. Droniou. Finite-volume schemes for noncoercive elliptic problems with Neumann boundary conditions. *IMA J. Numer. Anal.*, 31(1):61–85, 2011.

[14] C. Chainais-Hillairet and T. O. Gallouët. Study of a pseudo-stationary state for a corrosion model: existence and numerical approximation. *Nonlinear Anal. Real World Appl.*, 31:38–56, 2016.

[15] C. Chainais-Hillairet and I. Lacroix-Violet. On the existence of solutions for a drift-diffusion system arising in corrosion modelling. *DCDS-B*, 20(1):77–92, 2014.

[16] C. Chainais-Hillairet, B. Merlet, and A. Zurek. Convergence of a finite volume scheme for a parabolic system with a free boundary modeling concrete carbonation. *ESAIM Math. Model. Numer. Anal.*, 52(2):457–480, 2018.

[17] M. Chatard. Asymptotic behavior of the Scharfetter-Gummel scheme for the drift-diffusion model. In *Finite volumes for complex applications VI. Problems & perspectives. Volume 1, 2*, volume 4 of Springer Proc. Math., pages 235–243. Springer, Heidelberg, 2011.

[18] K. Deimling. *Nonlinear Functional Analysis*. Springer Berlin, Heidelberg, 1985.

[19] J. Donea, A. Huerta, J.-Ph. Ponthot, and A. Rodríguez-Ferran. Arbitrary Lagrangian–Eulerian methods. *Encyclopaedia of Computational Mechanics*, 2004.

[20] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. *The Gradient Discretisation Method*, volume 42 of *Mathématiques et Applications*. Springer International Publishing, Cham, 2018.

[21] 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.

[22] H. Gajewski and K. Gröger. Reaction-diffusion processes of electrically charged species. *Math. Nachr.*, 177:109–130, 1996.

[23] T. Gallouët and J.-C. Latché. Compactness of discrete approximate solutions to parabolic PDEs–Application to a turbulence model. *Comm. Pure Appl. Anal.*, 11(6), 2012.

[24] A. Glitzky. Analysis of electronic models for solar cells including energy resolved defect densities. *Math. Methods. App. Sci.*, 34(16):1980–1998, 2011.

[25] A. Glitzky. An electronic model for solar cells including active interfaces and energy resolved defect densities. *SIAM Journal on Mathematical Analysis*, 44(6):3874–3900, 2012.

[26] A. Glitzky, K. Gröger, and R. Hünlich. Free energy and dissipation rate for reaction diffusion processes of electrically charged species. *Appl. Anal.*, 60:201–217, 1996.
[27] A. M. Il'in. Differencing scheme for a differential equation with a small parameter affecting the highest derivative. *Math. Notes Acad. Sci. USSR*, 6:596-602, 1969.

[28] R. D. Lazarov, I. D. Mishev, and P. S. Vassilevski. Finite volume methods for convection-diffusion problems. *SIAM J. Numer. Anal.*, 33(1):31–55, 1996.

[29] B. Li, J. Lowengrub, A. Ratz, and A. Voigt. Geometric evolution laws for thin crystalline films: modeling and numerics. *Commun. Comput. Phys.*, 6(3):433, 2009.

[30] A. Mielke. A gradient structure for reaction-diffusion systems and for energy-drift-diffusion systems. *Nonlinearity*, 24(4):1329–1346, 2011.

[31] A. Mielke. Non-equilibrium steady states as saddle points and EDP-convergence for slow-fast gradient systems. *Journal of Mathematical Physics*, 64(12), 2023.

[32] A. Pimpinelli and J. Villain. *Physics of cristal growth*. Cambridge University Press, 1999.

[33] J. W. Portegies and M. A. Peletier. Well-posedness of a parabolic moving-boundary problem in the setting of Wasserstein gradient flows. *Interfaces Free Bound.*, 12(2):121–150, 2010.

[34] V. Rybalko and L. Berlyand. Emergence of traveling waves and their stability in a free boundary model of cell motility. *Trans. Amer. Math. Soc.*, 376(03):1799–1844, 2023.

[35] D.L. Scharfetter and H.K. Gummel. Large-signal analysis of a silicon Read diode oscillator. *IEEE Transactions on Electron Devices*, 16(1):64–77, 1969.