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

HAL
open science

# Study of a pseudo-stationary state for a corrosion model: existence and numerical approximation

Claire Chainais-Hillairet, Thomas Gallouët

► To cite this version:

Claire Chainais-Hillairet, Thomas Gallouët. Study of a pseudo-stationary state for a corrosion model: existence and numerical approximation. Nonlinear Analysis: Real World Applications, 2016, 31, pp.38-56. 10.1016/j.nonrwa.2016.01.010. hal-01147621

HAL Id: hal-01147621

https://inria.hal.science/hal-01147621

Submitted on 30 Apr 2015

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.

Public Domain
STUDY OF A PSEUDO-STATIONARY STATE FOR A
CORROSION MODEL: EXISTENCE AND NUMERICAL
APPROXIMATION

CLAIRE CHAINAIS-HILLAIRET AND THOMAS O. GALLOUËT

ABSTRACT. In this paper, we consider a system of partial differential equations describing the pseudo-stationary state of a dense oxide layer. We investigate the question of existence of a solution to the system and we design a numerical scheme for its approximation. Numerical experiments with real-life data shows the efficiency of the method. Then, the analysis is fulfilled on some simplified models.

# 1. INTRODUCTION

1.1. **General framework of the study.** The concept for long term storage of high-level radioactive waste in France under study is based on an underground repository. The waste shall be confined in a glass matrix and then placed into cylindrical steel canisters. These containers shall be placed into micro-tunnels in the highly impermeable Callovo-Oxfordian claystone layer at a depth of several hundred meters.

The long-term safety assessment of the geological repository has to take into account
the degradation of the carbon steel used for the waste overpacks and the cell dis-
posal liners, which are in contact with the claystone formation. This degradation is
mainly caused by generalized corrosion processes which take place under anaerobic
conditions. As a tool to investigate the corrosion processes at the surface of the
carbon steel canisters, the Diffusion Poisson Coupled Model (DPCM) for corrosion
has been developed by Bataillon et al. [4].

The DPCM model describes the evolution of a dense oxide layer in the region of contact between the claystones and the metal. This is based on a coupled system of electromigration – diffusion equations for the transport of the charge carriers in the oxide layer, and a Poisson equation for the electric potential. The interactions between the oxide layer and the clay soil or the metal are described in terms of Robin boundary conditions. The system includes moving boundary equations based on the Pilling-Bedworth ratio. As the oxide layer is very thin compared to the waste overpack size, it is sufficient to consider the model in one space dimension.

In [5], Bataillon *et al* proposed a numerical method for the approximation of the DPCM model. This numerical scheme has been implemented in the simulation code CALIPSO, developed at ANDRA (the french nuclear waste management agency). The long time simulations highlighted the existence of a pseudo stationary state:

<footnote>The authors acknowledge support from ANDRA, from Inria-Mephysto team, from ERC QUAN-
THOM and from Labex CEMPI (ANR-11-LABX-0007-01).</footnote>
a state where all the density profiles, as the profile of the electric potential and the size of the oxide layer, are constant, while both interfaces are moving at the same velocity. Moreover, it appears that the pseudo stationary state can be reached more or less quickly regarding the value of the pH. In many cases of application, the significant solution to the DPCM model may be the pseudo stationary one and in order to reduce the numerical costs, it is interesting to propose a scheme for its direct computation.

In this paper, we propose an efficient numerical method for the approximation of the pseudo stationary state of the DPCM model. However, we also investigate the question of existence of a pseudo steady state on some simplified models. The existence of a pseudo steady state is not so obvious and seems strongly related to the boundary conditions. Indeed, in [1, 2, 3], Aiki and Muntean consider a diffusion equation on a moving domain, where only one boundary is free. In this case, they prove that the thickness of the domain increases following a $\sqrt{t}$-law; there is no steady-state.

## 1.2. Presentation of the DPCM model.
The DPCM model was introduced in [4]. It consists of three drift-diffusion equations for the densities of charge carriers – ferric cations (P), electrons (N) and oxygen vacancies (C)–, coupled with one elliptic equation for the electric potential (Ψ) and two evolutive equations for the interfaces of the domain ($X_0, X_1$).

We consider here the dimensionless DPCM model. We will first present the equations and then give a sense to all the parameters involved in the equations. Let us already mention that the scaling in time is performed with respect to the characteristic time of the cations. The equations for the carrier densities P, N, C, as the boundary conditions, have the same form. For $U = P, N$ or $C$, they are written:

$$ (1a) \qquad \epsilon_U \partial_t U + \partial_x J_U = 0, \quad J_U = -\partial_x U - z_U U \partial_x \Psi \quad \text{in } (X_0(t), X_1(t)), \forall t \ge 0, $$

$$ (1b) \qquad -J_U - U X'_0(t) = r_U^0(U(X_0(t)), \Psi(X_0(t))) \quad \text{on } x = X_0(t), \forall t \ge 0, $$

$$ (1c) \qquad J_U + U X'_1(t) = r_U^1(U(X_1(t)), \Psi(X_1(t)), V) \quad \text{on } x = X_1(t), \forall t \ge 0, $$

where $z_U$ is the charge number of the carrier and $\epsilon_U$ is the ratio of the mobility coefficient with respect to the mobility coefficient of the cations (due to the choice of the time scaling). For $U = P, N, C$, we have respectively $z_U = 3, -1, 2$, and $\epsilon_U = 1, \frac{D_1}{D_2}, \frac{D_1}{D_3}$.

We also note that both functions $r_U^0$ and $r_U^1$ are linear and monotonically increasing with respect to their first argument. More precisely, the functions $r_U^0$ and $r_U^1$ have the following form:

$$ (2a) \qquad r_U^0(s, x) = \beta_U^0(x)s - \gamma_U^0(x), $$

$$ (2b) \qquad r_U^1(s, x, V) = \beta_U^1(V - x)s - \gamma_U^1(V - x), I $$
where $\beta_U^0, \beta_U^1, \gamma_U^0, \gamma_U^1$ are smooth positive functions defined as:

$$\gamma_P^0(x) = m_P^0 P^m e^{-3b_P^0 x}$$

$$\beta_P^0(x) = m_P^0 e^{-3b_P^0 x} + k_P^0 e^{3a_P^0 x}$$

$$\gamma_N^0(x) = m_N^0 e^{b_N^0 x}$$

$$\beta_N^0(x) = k_N^0 e^{-a_N^0 x}$$

$$\gamma_C^0(x) = m_C^0 e^{-2b_C^0 x}$$

$$\beta_C^0(x) = m_C^0 e^{-2b_C^0 x} + k_C^0 e^{2a_C^0 x}$$

$$\gamma_P^1(y) = k_P^1 P^m e^{3a_P^1 y}$$

$$\beta_P^1(y) = m_P^1 e^{-3b_P^1 y} + k_P^1 e^{3a_P^1 y}$$

$$\gamma_N^1(y) = k_N^1 N_{metal} \log(1 + e^{-y})$$

$$\beta_N^1(y) = m_N^1$$

$$\gamma_C^1(y) = k_C^1 e^{3a_C^1 y}$$

$$\beta_C^1(y) = m_C^1 e^{-3b_C^1 y} + k_C^1 e^{3a_C^1 y}$$

The equation on the electric potential $\Psi$ is a Poisson equation, written as:

$$(3a) \qquad -\lambda^2 \partial_{xx}^2 \Psi = 3P - N + 2C + \rho_{hl}, \quad x \in (X_0(t), X_1(t)),$$

$$(3b) \qquad \Psi - \alpha_0 \partial_x \Psi = \Delta \Psi_0^{pzc}, \quad x = X_0(t),$$

$$(3c) \qquad \Psi + \alpha_1 \partial_x \Psi = V - \Delta \Psi_1^{pzc}, \quad x = X_1(t).$$

The moving boundary equations are written as:

$$(4a) \qquad \frac{dX_0}{dt} = v_d^0(t) + \frac{dX_1}{dt} \left(1 - \frac{\Omega_{ox}}{m\Omega_{Fe}}\right),$$

$$(4b) \qquad \frac{dX_1}{dt} = -\frac{D_3}{4D_1} \frac{\Omega_{Fe}}{\Omega_{ox}} (J_C(X_1) + C X_1'(t)),$$

where $v_d^0(t)$ is the dissolution speed of the layer, given by $v_d^0(t) = k_d^0 e^{-5a_d^0 \Psi(X_0(t))}$.

The system is supplemented with initial conditions:

$$(5a) \qquad N(x, 0) = N^0(x), \quad P(x, 0) = P^0(x), \quad C(x, 0) = C^0(x), \quad x \in (0, 1),$$

$$(5b) \qquad X_0(0) = 0, \quad X_1(0) = 1.$$

Let us now shortly explain the parameters of the model. We first introduce the physical parameters:

* $D_1, D_2$ and $D_3$ are respectively the mobility or diffusion coefficients of cations, electrons and oxygen vacancies. $D_1$ and $D_3$ have the same order of magnitude, but $D_1 \ll D_2$ due to the difference of size between cations and electrons and the resulting difference of mobilities.

* $(a_U^0, b_U^0)$ for $U = P, N, C, r$, $(a_U^1, b_U^1)$ for $U = P, C$ and $a_d^0$ are positive transfer coefficients. They satisfy:

$$(6) \qquad a_U^0 + b_U^0 = 1 = a_U^1 + b_U^1 \quad \text{for } U = P, N, C.$$

* $\Omega_{ox}$ is the molar volume of the oxide, $\Omega_{Fe}$ is the molar volume of the metal and $m$ is the number of moles of iron per mole of oxide ($m = 3$ for magnetite).

The following parameters come from physical parameters but are concerned by the scaling. There are scaled parameters:
* $(m_P^i, k_P^i)_{i=0,1}, (m_N^i, k_N^i)_{i=0,1}, (m_C^i, k_C^i)_{i=0,1}, (n_N^0, p_N^0), k_d^0$ are interface kinetic functions. We assume that these functions are constant and strictly positive.

* $P^m$ is related to the maximum occupancy for octahedral cations in the host lattice.

* $N_{metal}$ is related to the electron density of state in the metal (Friedel model).

* $\rho_{hl}$ is the related to net charge density of the ionic species in the host lattice. We assume that $\rho_{hl}$ is homogeneous.

* $\Delta\Psi_0^{pzc}, \Delta\Psi_1^{pzc}$ are related respectively to the outer and the inner voltages of zero charge.

* $\lambda^2, \alpha_0, \alpha_1$ are positive dimensionless parameters coming from the scaling.

In the system (1)–(5), $V$ is an applied potential and the system corresponds to a “potentiostatic case”. However, the system (1)–(5) can be supplemented with an additional equation ensuring the electron charge balance at the inner interface :

$$ (7) \quad -3 \left( J_P + P X_1'(t) + \frac{D_3}{4D_1} (J_C + C X_1'(t)) \right) + \frac{D_2}{D_1} (J_N + N X_1'(t)) = \tilde{J}, \quad x = X_1(t) $$

This corresponds to a “galvanostatic case”. If $\tilde{J} = 0$, we speak of free corrosion and the resulting $V$ is called “free corrosion potential”.

## 1.3. Main results.
A numerical scheme was proposed in [5] in order to compute an approximate solution to the system (1)–(5) (or (1)–(7)). This scheme is a fully implicit Euler scheme in time and a finite volume scheme in space with a Scharfetter-Gummel approximation of the convection-diffusion fluxes. The numerical experiments, see for instance Figure 6.2 in [5] or Figure 36 in [4], show the apparition in large time of a pseudo steady state. The paper is devoted to the theoretical and numerical analysis of this steady state. It is structured as follows.

In Section 2, we define the steady state and we propose a finite volume scheme, closely related to the scheme from [5], in order to compute an approximate solution. We discuss the numerical implementation (as the scheme leads to a nonlinear system of equations). We also present and analyse some numerical results obtained with real-life data.

In Section 3, we investigate the question of existence of a steady state. We perform a complete analysis of two simplified models. The existence results are stated in Theorem 3.1 and Theorem 3.2. Then, Section 4 is devoted to the numerical analysis of the simplified models.

## 2. STUDY OF THE PSEUDO STATIONARY STATE

### 2.1. Definition.
The numerical experiments presented in [5, 4] show the apparition, in large time, of a pseudo stationary state. Indeed, it appears that the thickness of the oxide layer become constant, as well as the profiles of the densities and of the electric potential, while both interfaces are moving at the same velocity. In order to describe this steady state, we will use the following set of unknowns: the densities $P, N, C$, the electric potential $\Psi$, the thickness of the layer $\ell$ and the velocity of
the interfaces $\delta$. The set of equations will be deduced from the reformulation of the DPCM model on a fixed domain.

2.1.1. Reformulation of DPCM on fixed domain. In order to define the pseudo steady state, we reformulate the system (1)–(4) on a fixed domain, as in [5]. Therefore, we use the following change of variables:

$$ \bigcup_{0 \le t \le T} [X_0(t), X_1(t)] \times \{t\} \to [0, 1] \times [0, T] \\ (x, t) \mapsto \xi(x, t) = \left( \frac{x - X_0(t)}{X_1(t) - X_0(t)}, t \right) $$

We set

$$ (8) \qquad L(t) = X_1(t) - X_0(t). $$

The equations (1) on $U = P, N, C$ rewrite as:

$$ \varepsilon_U L(t) \partial_t (L(t) U) + \partial_\xi \hat{J}_U = 0, $$

$$ (9a) \qquad \hat{J}_U = -\partial_\xi U - [z_U \partial_\xi \Psi + \varepsilon_U L(t) (X_0'(t) + \xi L'(t))] U, \quad \xi \in (0, 1), $$

$$ (9b) \qquad -\hat{J}_U = L(t) r_U^0(U, \Psi), \quad \xi = 0, $$

$$ (9c) \qquad \hat{J}_U = L(t) r_U^1(U, \Psi, V), \quad \xi = 1. $$

The equations (3) on $\Psi$ rewrite as:

$$ (10a) \qquad -\frac{\lambda^2}{L(t)^2} \partial_{\xi\xi}^2 \Psi = 3P - N + 2C + \rho_{hl}, \quad \xi \in (0, 1), $$

$$ (10b) \qquad \Psi - \frac{\alpha_0}{L(t)} \partial_\xi \Psi = \Delta \Psi_0^{pzc}, \quad \xi = 0, $$

$$ (10c) \qquad \Psi + \frac{\alpha_1}{L(t)} \partial_\xi \Psi = V - \Delta \Psi_1^{pzc}, \quad \xi = 1. $$

Let us denote $\Pi = \frac{\Omega_{ox}}{m\Omega_{Fe}}$ and $K = \frac{D_3}{4D_1} \frac{\Omega_{Fe}}{\Omega_{ox}}$. Then, the moving boundary equations rewrite as:

$$ (11a) \qquad \frac{dX_0}{dt} = v_d^0(t) + \frac{dX_1}{dt}(1 - \Pi), $$

$$ (11b) \qquad \frac{dX_1}{dt} = -K r_C^1(C(1, t), \Psi(1, t), V), $$

$$ (11c) \qquad \text{with } v_d^0(t) = k_d^0 e^{-5a_d^0 \Psi(0, t)}. $$

Moreover, in the galvanostatic case, the additional equation becomes

$$ (12) \qquad -3 \left( \hat{J}_P + \frac{\varepsilon_C}{4} \hat{J}_C \right) + \varepsilon_N \hat{J}_N = \hat{J}, \quad \xi = 1. $$
2.1.2. The steady state model. We are now able to write the set of equations which will define the steady state model for DPCM. As it has already been written, the unknowns of this model are the densities of charge carriers P, N, C, the electric potential Ψ, the velocity of the interfaces δ and the thickness of the oxide layer ℓ. The set of equations is obtained from (9)–(11) by letting down the dependency with respect to time and by setting L(t) = ℓ and X'_0(t) = X'_1(t) = δ. It writes as:

• Equations and boundary conditions for U = P, N, C :

$$
(13a) \quad \partial_{\xi} \hat{J}_{U} = 0, \quad \hat{J}_{U} = -\partial_{\xi} U - (z_{U} \partial_{\xi} \Psi + \varepsilon_{U} \ell \delta) U, \quad \xi \in (0, 1),
$$

$$
(13b) \quad -\hat{J}_U = \ell r_U^0 (U, \Psi), \quad \xi = 0,
$$

$$
(13c) \quad \hat{J}_U = \ell r_U^1 (U, \Psi, V), \quad \xi = 1.
$$

• Equation and boundary condition for Ψ :

$$
(14a) \quad -\frac{\lambda^2}{\ell^2} \partial_{\xi\xi}^2 \Psi = 3P - N + 2C + \rho_{hl}, \quad \xi \in (0, 1),
$$

$$
(14b) \quad \Psi - \frac{\alpha_0}{\ell} \partial_{\xi} \Psi = \Delta \Psi_0^{pzc}, \quad \xi = 0,
$$

$$
(14c) \quad \Psi + \frac{\alpha_1}{\ell} \partial_{\xi} \Psi = V - \Delta \Psi_1^{pzc}, \quad \xi = 1.
$$

• Equation for the velocity δ and the thickness ℓ :

$$
(15a) \qquad \delta = \frac{v_d^0}{\Pi},
$$

$$
(15b) \qquad \ell = -K \frac{\hat{J}_C(1)}{\delta}
$$

$$
(15c) \quad \text{with } v_d^0 = k_d^0 e^{-5a_d^0 \Psi(0)}
$$

The set of equations (13)–(15) describing the steady state will be denoted (M) in
the sequel. It describes the potentiostatic case. In the galvanostatic case, V is
defined by an additional equation:

$$
(16) \qquad -3 \left( \hat{J}_P + \frac{\varepsilon_C}{4} \hat{J}_C \right) + \varepsilon_N \hat{J}_N = 0.
$$

The system (13)–(16) will be denoted (Mg).

2.2. **Numerical scheme.** In order to compute an approximate solution to (M) or (Mg), we propose a finite volume scheme close to the one introduced in [5]. Therefore, we consider a mesh for the domain [0, 1], which is not necessarily uniform, i.e a family of given points $(x_i)_{0\le i\le I+1}$ satisfying

$$
x_0 = 0 < x_1 < x_2 < \dots < x_I < x_{I+1} = 1.
$$

Then, for $1 \le i \le I - 1$, we define $x_{i+\frac{1}{2}} = \frac{x_i + x_{i+1}}{2}$ and we set $x_{\frac{1}{2}} = x_0 = 0$, $x_{I+\frac{1}{2}} = x_{I+1} = 1$. The cells of the mesh are the intervals $(x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}})$ for $1 \le i \le I$.
Let us set

$$h_i = x_{i+\frac{1}{2}} - x_{i-\frac{1}{2}}, \quad \text{for } 1 \le i \le I,$$

$$h_{i+\frac{1}{2}} = x_{i+1} - x_i, \quad \text{for } 0 \le i \le I$$

and $h = \max\{h_i, 1 \le i \le I\}$ is the size of the mesh.

The unknowns of the scheme for (M) are the densities $(N_i, P_i, C_i)_{0 \le i \le I+1}$ and the electric potential $(\Psi_i)_{0 \le i \le I+1}$, the velocity of the interfaces $\delta^h$ and the thickness of the domain $\ell^h$.

The scheme writes as:

• Scheme for $\Psi$:

$$(17a) \qquad -\frac{\lambda^2}{\ell^{h2}}(d\Psi_{i+\frac{1}{2}} - d\Psi_{i-\frac{1}{2}}) = h_i(3P_i - N_i + 2C_i + \rho_{hl}), \quad 1 \le i \le I,$$

$$(17b) \qquad \text{with } d\Psi_{i+\frac{1}{2}} = \frac{\Psi_{i+1} - \Psi_i}{h_{i+\frac{1}{2}}}, \quad 0 \le i \le I,$$

$$(17c) \qquad \Psi_0 - \frac{\alpha_0}{\ell^h} d\Psi_{\frac{1}{2}} = \Delta\Psi_0^{pzc},$$

$$(17d) \qquad \Psi_{I+1} + \frac{\alpha_1}{\ell^h} d\Psi_{I+\frac{1}{2}} = V - \Delta\Psi_1^{pzc}.$$

• Scheme for $U = P, N, C$:

$$(18a) \qquad \mathcal{G}_{U,i+\frac{1}{2}} - \mathcal{G}_{U,i-\frac{1}{2}} = 0, \quad 0 \le i \le I,$$

$$(18b) \qquad \mathcal{G}_{U,i+\frac{1}{2}} = \frac{1}{h_{i+\frac{1}{2}}} \left( B \left( h_{i+\frac{1}{2}} \left( z_U d\Psi_{i+\frac{1}{2}} + \varepsilon_U \ell^h \delta^h \right) \right) U_i - \right.$$

$$(18c) \qquad B \left( -h_{i+\frac{1}{2}} \left( z_U d\Psi_{i+\frac{1}{2}} + \varepsilon_U \ell^h \delta^h \right) \right) U_{i+1} \right), \quad 1 \le i \le I,$$

$$(18d) \qquad -\mathcal{G}_{U,\frac{1}{2}} = \ell^h r_U^0 (U_0, \Psi_0),$$

$$(18e) \qquad \mathcal{G}_{U,I+\frac{1}{2}} = \ell^h r_U^1 (U_{I+1}, \Psi_{I+1}),$$

where $B$ is given by the Bernoulli function

$$(19) \qquad B(x) = \frac{x}{e^x - 1} \quad \forall x \neq 0, \quad B(0) = 1.$$

• Scheme for $\delta$ and $\ell$:

$$(20a) \qquad \delta^h = \frac{k_d^0}{\Pi} e^{-5a_d^0 \Psi_0},$$

$$(20b) \qquad \ell^h = -K \frac{\mathcal{G}_{C,I+\frac{1}{2}}}{\delta^h}.$$

The scheme (17)–(20) will be denoted by (S). In the galvanostatic case, we have to discretize the additional equation (16). We write:

$$(21) \qquad -3 \left( \mathcal{G}_{P,I+\frac{1}{2}} + \frac{\varepsilon_C}{4} \mathcal{G}_{C,I+\frac{1}{2}} \right) + \varepsilon_N \mathcal{G}_{N,I+\frac{1}{2}} = 0.$$
It adds one unknown and one nonlinear equation to the previous system of nonlinear equations (S). In the sequel, we will denote this scheme (Sg).

Let us remark that the choice of the Bernoulli function for B corresponds to a Scharfetter-Gummel approximation of the convection-diffusion fluxes. These numerical fluxes have been introduced by Il'in in [12] and Scharfetter and Gummel in [15] for the numerical approximation of the drift-diffusion system arising in semiconductor modelling. It has been established by Lazarov, Mishev and Vassilevsky in [13] that they are second-order accurate in space. Dissipativity of the Scharfetter-Gummel scheme with a backward Euler time discretization for the classical drift-diffusion system was proved in [11] and Chatard in [10]. One crucial property of the Scharfetter-Gummel fluxes is that they generally preserve steady-states. We will see in Section 4 that this remains true in our case.

## 2.3. Implementation and numerical results.

### 2.3.1. Implementation.
The scheme (S) can be seen as the limit $\Delta t \to +\infty$ of the scheme proposed in [5] (in this paper, the system is the evolutive DPCM model and $\Delta t$ is the time step of the scheme). Therefore, the implementation is similar to the one proposed in [5]. The main difference is that the equations for the moving boundary interfaces, written on $X_0$ and $X_1$ in the evolutive case, become equations on the velocity $\delta^h$ and the thickness $\ell^h$.

The scheme (S) leads to a nonlinear system of equations $F(Z) = 0$, whose unknown is

$$Z = ((\Psi_i, P_i, N_i, C_i)_{0\le i\le I+1}, \delta^h, \ell^h) \in \mathbb{R}^{4(I+2)+2}$$

in the potentiostatic case or $Z = ((\Psi_i, P_i, N_i, C_i)_{0\le i\le I+1}, \delta^h, \ell^h, V) \in \mathbb{R}^{4(I+2)+3}$ in the galvanostatic case. A Newton's method is used to solve the nonlinear system of equations, with a Schur complement technique for the solution of linear systems. We refer to [5] for more details. The implementation of the scheme (S) has been done in the code CALIPSO and we will denote from now on CALIPSO-(S) this part of the code.

One difficult task is the choice of the initialization of the Newton's method. We can start with a set of values for $((P_i), (N_i), (C_i))_{0\le i\le I}$ given in the input data file and compute $(\Psi_i)_{0\le i\le I+1}$ as a solution to (17). In this case, we choose some constant values for the densities and the convergence of the Newton's method is not always ensured. In practice, it works only for very few values of the pH. With the set of parameters given in Appendix A, it works for instance for pH=9.3.

Then, an other strategy is to apply a continuation process: we start with a pH value that is supporting the first method and we use the computed solution as initial data for a new value of pH, close to the preceding one. Thus, using small increments for the value of the pH, we are able to compute the pseudo stationary state for any value of the pH, in the potentiostatic as in the galvanostatic case.

### 2.3.2. Numerical experiments.

#### Profiles.
We first want to compare the different profiles computed with the code CALIPSO for the evolutive system at a large time and the pseudo steady state directly computed with the code CALIPSO-(S). In practice, we consider that the steady state is reached by the evolutive scheme when the relative variation of the thickness of the oxide layer is small with respect to a given parameter named $\epsilon_{tol}$. This parameter must be sufficiently small in order to obtain accurate results. In practice, we choose $\epsilon_{tol} = 10^{-6}$.

The set of parameters defining the test case is given in Appendix A. We choose a Tchebychev mesh (so that the mesh is refined at the boundaries, see [5]) with 2000 cells. In the evolutive case, the time step is $\Delta t = 1000$s (before scaling).

Figures 2.1 and 2.3 show the profiles of the densities and of the electric potential in the potentiostatic case with $V = 0.5$ and pH=9.3. We observe that the profiles computed with the evolutive scheme and with (S) are similar. The conclusion is the same in the galvanostatic case with pH=9.3, see Figures 2.3 and 2.4.

<figure><img src="image_2.png" /><figcaption>FIGURE 2.1. Potentiostatic case with $V = 0.5$, at pH=9.3. Density profiles at the steady-state, (a) reached by the code CALIPSO for the evolutive system, (b) directly computed by CALIPSO-(S).</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIGURE 2.2. Potentiostatic case with V = 0.5, at pH=9.3. Profile of the electric potential at the steady-state, (a) reached by the code CALIPSO for the evolutive system, (b) directly computed by CALIPSO-(S).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2.3. Galvanostatic case at pH=9.3. Density profiles at the steady-state, (a) reached by the code CALIPSO for the evolutive system, (b) directly computed by CALIPSO-(S).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIGURE 2.4. Galvanostatic case at pH=9.3. Profile of the electric potential at the steady-state, (a) reached by the code CALIPSO for the evolutive system, (b) directly computed by CALIPSO-(S).</figcaption></figure>

## Efficiency of the direct computation with (S) or (Sg)

We now compare the codes CALIPSO and CALIPSO-(S) in terms of the number of iterations of the Newton's method. Indeed, as the scheme for the evolutive DPCM model implemented in the code CALIPSO is a fully implicit scheme, the discrete solution at each time step is computed with a Newton's method. In this case, the values for the initialization of the method are given by the solution at the preceding time step.

In Table 2.1, we present the number of iterations of the Newton's method necessary for the computation of the steady state with CALIPSO and CALIPSO-(S) for two values of pH. At pH=9.3, we can initialize the Newton's method in CALIPSO-(S) with a constant set of discrete densities. But, the solution at pH=10.2 is obtained thanks to the continuation process described above: starting from pH=9.3 and increasing the value of pH up to 10.2 by increments of pH equal to $10^{-2}$. The stopping criterion in both codes are the same and equal to $\epsilon_{Newt} = 10^{-8}$. Table 2.1 shows the efficiency of the direct computation of the steady state with (S) or (Sg). Figure 2.5 shows the dependency of the thickness of the oxide layer with respect to the pH value. We remark that the pseudo steady state is reached earlier for low values of the pH (close to 9) than for high values. Therefore, the computation of the steady state with CALIPSO-(S) is the more useful than the pH is high.

<table><caption>TABLE 2.1. Number of Newton iterations</caption><thead><tr><th rowspan="2"></th><th colspan="2">Potentiostatic case</th><th colspan="2">Galvanostatic case</th></tr><tr><th>pH = 9.3</th><th>pH = 10.2</th><th>pH = 9.3</th><th>pH = 10.2</th></tr></thead><tbody><tr><th>CALIPSO-(S)</th><td>12</td><td>57</td><td>23</td><td>68</td></tr><tr><th>CALIPSO</th><td>2359</td><td>12037</td><td>97</td><td>3654</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>FIGURE 2.5. Dependency of the thickness of the oxide layer with respect to the pH value. Reached by CALIPSO-(S) and CALIPSO for increasing final time, (a) potentiostatic case, (b) galvanostatic case.</figcaption></figure>

Sensitivity of the models (M) and (Mg) to the parameters.

The thickness of the oxide layer at the steady state strongly depends on the pH and, in the potentiostatic case, on the applied potential V. We can already see on Figure 2.5 that the steady state thickness follows a logarithmic law with respect to pH with a slope equal to 1. It is related to the fact that $k_d^0$ is proportional to $10^{-pH}$. The depassivation pH is defined as the supremum of the pH values for which the thickness of the corrosion layer is smaller the 1nm.

Figure 2.6 shows that the depassivation pH highly depends on the value of $k_d^0$ and therefore the importance of giving a very good approximation of this data. The slope of the logarithmic law does not change. The value of $k_{d,ref}^0$ is given in the appendix A

### 3. MATHEMATICAL STUDY OF SIMPLIFIED MODELS

The existence of a solution to the evolutive DPCM model (1)–(5) is an open question. The main difficulties are of three types: the system is strongly coupled (the coupling arises in the equations and in the boundary conditions), the boundary conditions are Robin boundary conditions and the interfaces are moving. In [9], Chainais-Hillairet and Lacroix-Violet proves the existence of a solution in a simplified case where only electrons and cations are taken into account and therefore the domain is fixed. Convergence of backward Euler scheme in time and finite volume scheme in space for the same simplified system has been established in [7]. Chainais-Hillairet and Lacroix-Violet have also proved in [8] the existence of a steady state for the two species system on a fixed domain, while the convergence of a finite volume scheme is studied in [6].
<figure><img src="image_7.png" /><figcaption>FIGURE 2.6. Dependency of the thickness of the oxide layer with respect to the pH value. Reached by CALISPO-(S) for multiple values of $k_d^0$. (a) potentiostatic case, (b) galvanostatic case.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 2.7. Dependency of the thickness of the oxide layer with respect to the applied potential V. Reached by CALISPO-(S) for multiple values of the pH in the potentiostatic case.</figcaption></figure>

In this paper, we have introduced the system defining the steady state of the full DPCM model: (M) in the potentiostatic case and (Mg) in the galvanostatic case. From now on, we will focus on the potentiostatic case. The coupling is still strong. However, the numerical experiments done for (M) show, for different values of the pH, that the electric potential Ψ has a linear profile, see Figure 3.1. This observation has led us to some simplifications in (M), in order to get first theoretical results.

We will consider two simplified models: (TM1) and (TM2). The model (TM1) is obtained by setting $Ψ = 0$ in (M). This assumption on Ψ has no physical sense (see Figure 3.1). However, the study of (TM1) will permit us to highlight some hypotheses on the data needed for the existence of a solution.

The simplified model (TM2) is defined by cancelling the right hand side of the equation on Ψ (14). Therefore Ψ is effectively an affine function on [0, 1] which is explicitely determined, but the densities and the electric potential are still coupled via the boundary conditions in the equations for the densities.
<figure><img src="image_9.png" /><figcaption>FIGURE 3.1. Profiles of the electric potential for different values of pH in the potentiostatic case.</figcaption></figure>

In this Section, we first introduce the simplified models (TM1) and (TM2). Then, we prove the existence and uniqueness of a solution for both models. The question of the positivity of the densities and of the thickness of the oxide layer is also investigated. The numerical analysis of the schemes deduced from (S) for both simplified models is proposed in Section 4.

### 3.1. Presentation of the simplified models.

*First simplified model: (TM1).* We first focus on the boundary issue, by fixing $\Psi = 0$. The model (TM1) is given by the following set of equations:

* Equation and boundary conditions for C:

$$ (22a) \qquad \partial_{\xi} \hat{J}_{C} = 0, \quad \hat{J}_{C} = -\partial_{\xi} C - \epsilon_{C} \ell \delta C, \quad \xi \in (0, 1), $$

$$ (22b) \qquad -\hat{J}_{C} = \ell r_{C}^{0}(C, 0), \quad \xi = 0, $$

$$ (22c) \qquad \hat{J}_{C} = \ell r_{C}^{1}(C, 0, V), \quad \xi = 1. $$

* Equation for the velocity $\delta$ and the thickness $\ell$:

$$ (23a) \qquad \delta = \frac{k_{d}^{0}}{\Pi}, $$

$$ (23b) \qquad \ell = -K \frac{\hat{J}_{C}(1)}{\delta}. $$

* Equation and boundary conditions for $U=P, N$.

$$ (24a) \qquad \partial_{\xi} \hat{J}_{U} = 0, \quad \hat{J}_{U} = -\partial_{\xi} U - \epsilon_{U} \ell \delta U, \quad \xi \in (0, 1), $$

$$ (24b) \qquad -\hat{J}_{U} = \ell r_{U}^{0}(U, 0), \quad \xi = 0, $$

$$ (24c) \qquad \hat{J}_{U} = \ell r_{U}^{1}(U, 0, V), \quad \xi = 1. $$

Let us just remark that (24) is in practice independent from (22)-(23). Indeed, if we prove existence and uniqueness of a solution to (22)-(23), we will directly get existence and uniqueness of a solution to (22)-(24). Therefore, the model (TM1) corresponds to (22)-(23).
Second simplified model: (TM2). The simplified model (TM2) is obtained by can-
celling the right hand side of the equation on Ψ (14). It is given by the following
set of equations:

• Equation and boundary conditions for Ψ:

$$
(25a) \qquad -\frac{\lambda^2}{\ell^2} \partial_{xx}^2 \Psi = 0, \quad x \in (X_0(t), X_1(t)),
$$

$$
(25b) \qquad \Psi - \frac{\alpha_0}{\ell} \partial_x \Psi = \Delta \Psi_0^{pzc}, \quad x = X_0(t),
$$

$$
(25c) \qquad \Psi + \frac{\alpha_1}{\ell} \partial_x \Psi = V - \Delta \Psi_1^{pzc}, \quad x = X_1(t).
$$

• Equation and boundary conditions for C

$$
(26a) \quad \partial_{\xi} \hat{J}_{C} = 0, \quad \hat{J}_{C} = -\partial_{\xi} C - (z_{C} \partial_{\xi} \Psi + \epsilon_{C} \ell \delta C), \quad \xi \in (0, 1),
$$

$$
(26b) \quad -\hat{J}_C = \ell r_C^0(C, \Psi(0)), \quad \xi = 0,
$$

$$
(26c) \quad \hat{J}_C = \ell r_C^1(C, \Psi(1), V), \quad \xi = 1.
$$

• Equation for the velocity $\delta$ and the thickness $\ell$:

$$
(27a) \qquad \delta = \frac{k_d^0 e^{-5a_d^0 \Psi(0)}}{\Pi},
$$

$$
(27b) \qquad \ell = -K \frac{\hat{J}_C(1)}{\delta}.
$$

• Equation and boundary conditions for U = P, N.

$$
(28a) \quad \partial_{\xi} \hat{J}_{U} = 0, \quad \hat{J}_{U} = -\partial_{\xi} U - (z_{U} \partial_{\xi} \Psi + \varepsilon_{U} \ell \delta U), \quad \xi \in (0, 1),
$$

$$
(28b) \quad -\hat{J}_U = \ell r_U^0 (U, \Psi(0)), \quad \xi = 0,
$$

$$
(28c) \quad \hat{J}_U = \ell r_U^1 (U, \Psi(1), V), \quad \xi = 1.
$$

As for the model (TM1), the computation of N and P, solutions to (28), is decou-
pled from (25)–(27). Therefore, the model (TM2) corresponds to (25)–(27).

3.2. Existence and uniqueness of a solution for (TM1). Let us first remark that (22a) implies that $\hat{J}_C$ is constant. Therefore, (22) can be rewritten as the following boundary value problem in (0, 1):

$$
(29a) \qquad \partial_{\xi} C + \varepsilon_C \ell \delta C = -\hat{J}_C, \quad \xi \in (0, 1)
$$

$$
(29b) \qquad C(0) = \frac{1}{\beta_C^0(0)} \left( \gamma_C^0(0) - \frac{\hat{J}_C}{\ell} \right)
$$

$$
(29c) \qquad C(1) = \frac{1}{\beta_C^1(V)} \left( \gamma_C^1(V) + \frac{\hat{J}_C}{\ell} \right).
$$
By standard arguments for linear differential equations, we find C solution to (29a)-(29b) and then deduce $\hat{J}_C$ from (29c). It yields:

$$ (30a) \quad C(\xi) = \left( C(0) + \frac{\hat{J}_C}{\varepsilon_C \ell \delta} \right) e^{-\varepsilon_C \ell \delta \xi} - \frac{\hat{J}_C}{\varepsilon_C \ell \delta}, $$

$$ (30b) \quad -\hat{J}_C = \left( \frac{\gamma_C^1(V)}{\beta_C^1(V)} - \frac{\gamma_C^0(0)}{\beta_C^0(0)} e^{-\varepsilon_C \ell \delta} \right) \left[ \frac{1}{\ell \beta_C^1(V)} + \frac{e^{-\varepsilon_C \ell \delta}}{\ell \beta_C^0(0)} + \frac{1 - e^{-\varepsilon_C \ell \delta}}{\varepsilon_C \ell \delta} \right]^{-1}. $$

The velocity of the boundary interfaces $\delta$ is explicitly given by (23a). The difficulty is to find the value of $\ell$ satisfying (23a). Let us remark that (23) implies:

$$ \frac{\hat{J}_C}{\ell \delta} = -\frac{1}{K} \quad \text{and} \quad \frac{\hat{J}_C}{\ell} = -\frac{k_d^0}{\Pi K}. $$

Plugging the first equality into the boundary conditions (29b) and (29c), we get:

$$ (31) \quad \begin{cases} C(0) = \frac{1}{\beta_C^0(0)} & \left( \gamma_C^0(0) + \frac{k_d^0}{\Pi K} \right), \\ C(1) = \frac{1}{\beta_C^1(V)} & \left( \gamma_C^1(V) - \frac{k_d^0}{\Pi K} \right), \end{cases} $$

Then, writing (30a) at $\xi = 1$ and using the second equality yields:

$$ (32) \quad \ell = \frac{\Pi}{\varepsilon_C k_d^0} \ln \left( \frac{\varepsilon_C K C(0) - 1}{\varepsilon_C K C(1) - 1} \right). $$

We have now to verify if $\ell$, given by (32), is well defined and to check if $\ell > 0$. Thanks to (32), the existence of $\ell > 0$ is ensured if and only if one of these two cases appears :

$$ (i) \text{ either } [\varepsilon_C K C(0) - 1 > \varepsilon_C K C(1) - 1 > 0] $$

$$ (ii) \text{ or } [\varepsilon_C K C(0) - 1 < \varepsilon_C K C(1) - 1 < 0]. $$

Condition (i) reduces to

$$ (33) \quad \varepsilon_C K C(1) > 1 \quad \text{and} \quad C(0) > C(1). $$

It leads to

$$ \varepsilon_C K \gamma_C^1(V) - \beta_C^1(V) > \varepsilon_C \frac{k_d^0}{\Pi} \quad \text{and} \quad \beta_C^1(V) \left( \gamma_C^0(0) + \frac{k_d^0}{\Pi K} \right) \geq \beta_C^0(0) \left( \gamma_C^1(V) - \frac{k_d^0}{\Pi K} \right), $$

or equivalently

$$ K \gamma_C^1(V) - \frac{\beta_C^1(V)}{\varepsilon_C} > \frac{k_d^0}{\Pi} > K \frac{\beta_C^0(0) \gamma_C^1(V) - \beta_C^1(V) \gamma_C^0(0)}{\beta_C^0(0) + \beta_C^1(V)}. $$

Similarly, condition (ii) reduces to

$$ K \gamma_C^1(V) - \frac{\beta_C^1(V)}{\varepsilon_C} < \frac{k_d^0}{\Pi} < K \frac{\beta_C^0(V) \gamma_C^1(V) - \beta_C^1(V) \gamma_C^0(V)}{\beta_C^0(V) + \beta_C^1(V)}. $$
Setting

$$ (34a) \qquad K_1 = \Pi \left( K \gamma_C^1(V) - \frac{\beta_C^1(V)}{\varepsilon_C} \right) $$

$$ (34b) \qquad K_2 = \Pi K \frac{\beta_C^0(0)\gamma_C^1(V) - \beta_C^1(V)\gamma_C^0(0)}{\beta_C^0(0) + \beta_C^1(V)}, $$

we find that the existence of $\ell > 0$ is ensured if and only if

$$ (35) \qquad \min(K_1, K_2) < k_d^0 < \max(K_1, K_2). $$

**Theorem 3.1.** Let assume $[\min(K_1, K_2) < k_d^0 < \max(K_1, K_2)]$, where $K_1$ and $K_2$ are defined by (34). Then, the problem (TM1) has a unique solution $(C_s, \ell_s, \delta_s, ) \in C^\infty([0, 1], \mathbb{R}) \times (\mathbb{R}_+^*)^2$ defined by:

$$ (36a) \qquad \delta_s = \frac{k_d^0}{\Pi}, $$

$$ (36b) \qquad \ell_s = \frac{\Pi}{\varepsilon_C k_d^0} \ln \left( \frac{\varepsilon_C K C_0 - 1}{\varepsilon_C K C_1 - 1} \right), $$

$$ (36c) \qquad C_s(\xi) = \left( C_0 - \frac{1}{\varepsilon_C K} \right) e^{-\varepsilon_C \ell_s \delta_s \xi} + \frac{1}{\varepsilon_C K}, $$

with $C_0 = \frac{1}{\beta_C^0(0)} \left( \gamma_C^0(0) + \frac{k_d^0}{\Pi K} \right)$ and $C_1 = \frac{1}{\beta_C^1(V)} \left( \gamma_C^1(V) - \frac{k_d^0}{\Pi K} \right)$.

Moreover

$$ (37) \qquad 0 \le \min \left( C_0, \frac{1}{\varepsilon_C K} \right) \le C_s(\xi) \le \max \left( C_0, \frac{1}{\varepsilon_C K} \right) \quad \forall \xi \in (0, 1). $$

*Proof.* We have already proved the existence and uniqueness of the solution, with $\ell_s > 0$. It remains to prove the bounds for $C_s$. We note that $C_s(0) = C_0$ is positive and that $C_s$ defined by (36c) is monotone. More precisely, in the case (i) described above, $C_s$ is decreasing and $1/(\varepsilon_C K) < C_1 \le C_s(\xi) \le C_0$ for all $\xi \in [0, 1]$, while in the case (ii), $C_s$ is increasing and $C_0 \le C_s(\xi) \le C_1 \le 1/(\varepsilon_C K)$ for all $\xi \in [0, 1]$. Finally, (37) holds in both cases, which concludes the proof.

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

**Remark 3.1.** When the solution to (TM1) is known, it is clear that the system (24) has a unique solution $U_s$ for $U = P$ or $N$. This solution is given by

$$ (38) \qquad \begin{cases} -\hat{J}_U = \left( \frac{\gamma_U^1(V)}{\ell_s \beta_U^1(V)} - \frac{\gamma_U^0(0)}{\ell_s \beta_U^0(0)} e^{-\varepsilon_U \ell_s \delta_s} \right) \left[ \frac{1}{\ell_s \beta_U^1(V)} + \frac{e^{-\varepsilon_U \ell_s \delta_s}}{\ell_s \beta_U^0(0)} + \frac{1 - e^{-\varepsilon_U \ell_s \delta_s}}{\varepsilon_U \ell_s \delta_s} \right]^{-1}, \\ U_s(\xi) = \left( \frac{\gamma_U^0(0)}{\beta_U^0(0)} - \frac{\hat{J}_U}{\ell_s \beta_U^0(0)} + \frac{\hat{J}_U}{\varepsilon_U \ell_s \delta_s} \right) e^{-\varepsilon_U \ell_s \delta_s \xi} - \frac{\hat{J}_U}{\varepsilon_U \ell_s \delta_s} \quad \forall \xi \in [0, 1]. \end{cases} $$

It is possible to impose new conditions on the data in order to ensure that $\hat{J}_U$ are nonpositive or that $P_s$ and $N_s$ are nonnegative.
**Remark 3.2.** The parameter $k_d^0$ is necessarily positive. Therefore, Theorem 3.1 is worthy of interest only if $\max(K_1, K_2) > 0$. This condition is satisfied either if $\beta_C^0(0)\gamma_C^1(V) \ge \beta_C^1(V)\gamma_C^0(0)$ or if $K\varepsilon_C\gamma_C^1(V) \ge \beta_C^1(V)$. Using the definition of the functions $\beta_C^0, \beta_C^1, \gamma_C^0, \gamma_C^1$ given in Section 1.2 and the hypothesis (6), we can see that the condition is satisfied if

$$e^{3V} \ge \min \left( \frac{m_C^0 m_C^1}{k_C^0 k_C^1}, \frac{m_C^1}{k_C^1} \frac{1}{\varepsilon_C K - 1} \right)$$

3.3. **Existence of a solution for (TM2).** In order to prove the existence of a solution for (TM2), we start with a given $\ell > 0$ and study the system (25), (26), (27a), see Lemma 3.1. Afterwards, we will prove the existence of $\ell$ satisfying (27b) thanks to an intermediate value theorem.

**Lemma 3.1.** Let $\ell > 0$. There exists a unique solution $(\Psi, C, \delta) \in (C^\infty[0, 1])^2 \times \mathbb{R}_+^*$ to the system (25)-(26)-(27a).

*Proof.* In the model (TM2), the electric potential is no more equal to 0. However, for a given $\ell > 0$, we can still compute $\Psi$ solution to (25). We have:

$$ \begin{array}{l} \Psi(\xi) = (\Psi(1) - \Psi(0))\xi + \Psi(0), \\ (39) \qquad \Psi(0) = \frac{(\ell + \alpha_1) \Delta\Psi_0^{pzc} + \alpha_0 (V - \Delta\Psi_1^{pzc})}{\ell + \alpha_0 + \alpha_1}, \\ \qquad \Psi(1) = \frac{\alpha_1 \Delta\Psi_0^{pzc} + (\ell + \alpha_0) (V - \Delta\Psi_1^{pzc})}{\ell + \alpha_0 + \alpha_1}. \end{array} $$

As $\max(\Delta\Psi_0^{pzc}, V - \Delta\Psi_1^{pzc})$ and $\min(\Delta\Psi_0^{pzc}, V - \Delta\Psi_1^{pzc})$ are sur- and sub-solutions of (25), $\Psi$ is bounded independently of $\ell$:

$$(40) \quad \min(\Delta\Psi_0^{pzc}, V - \Delta\Psi_1^{pzc}) \le \Psi(\xi) \le \max(\Delta\Psi_0^{pzc}, V - \Delta\Psi_1^{pzc}), \quad \forall \xi \in (0, 1).$$

According to formula (39), we remark that the solution to (25) tends to a constant function when $\ell$ tends to 0, where the constant is

$$\tilde{\Psi} = \frac{\alpha_1 \Delta\Psi_0^{pzc} + \alpha_0 (V - \Delta\Psi_1^{pzc})}{\alpha_0 + \alpha_1}.$$

When $\Psi$ is defined, we can deduce $\delta$:

$$(41) \qquad \delta = \frac{k_d^0}{\Pi} e^{-5a_d^0 \Psi(0)} > 0.$$

Then, knowing $\ell$, $\Psi$ and $\delta$, it is now possible to compute $C$ solution to (26) (and similarly $P$ and $N$ solutions to (28)). Therefore, we use the Slotboom's change of variables, which is classical for drift-diffusion equations, see for instance [14]. Defining $u(\xi) = U(\xi)e^{z_U\Psi(\xi)+\varepsilon_U\ell\delta\xi}$ for $U = C$ or $U = P, N$ we can rewrite (26) and (28) as:

$$(42a) \qquad \partial_{\xi} \hat{J}_U = 0, \quad \hat{J}_U = -e^{-z_U \Psi(\xi) - \varepsilon_U \ell \delta \xi} \partial_{\xi} u, \quad \xi \in (0, 1),$$

$$(42b) \qquad \partial_{\xi} u - \ell \beta_U^0 (\Psi(0)) u = -\ell \gamma_U^0 (\Psi(0)) e^{z_U \Psi(0)}, \quad \xi = 0,$$

$$(42c) \qquad \partial_{\xi} u + \ell \beta_U^1 (V - \Psi(1)) u = \ell \gamma_U^1 (V - \Psi(1)) e^{z_U \Psi(1) + \varepsilon_U \ell \delta}, \quad \xi = 1,$$
Let us set $B_U^0 = \beta_U^0(\Psi(0))$, $B_U^1 = \beta_U^1(V - \Psi(1))$, $\Gamma_U^0 = \gamma_U^0(\Psi(0))$ and $\Gamma_U^1 = \gamma_U^1(V - \Psi(1))$. Following the lines of the proof of Proposition 1 in [8], we obtain:

(43a) 	$u(\xi) = u(0) - \hat{J}_U \int_0^\xi e^{z_U \Psi(\zeta) + \varepsilon_U \ell \delta \zeta} d\zeta$

(43b) 	$-\hat{J}_U = \frac{\Gamma_U^1 e^{z_U \Psi(1) + \varepsilon_U \ell \delta} / B_U^1 - \Gamma_U^0 e^{z_U \Psi(0)} / B_U^0}{e^{z_U \Psi(0)} / (\ell B_U^0) + e^{z_U \Psi(1) + \varepsilon_U \ell \delta} / (\ell B_U^1) + \int_0^1 e^{z_U \Psi(\zeta) + \varepsilon_U \ell \delta \zeta} d\zeta}$

(43c) 	$u(0) = \Gamma_U^0 e^{z_U \Psi(0)} / B_U^0 - \hat{J}_U e^{z_U \Psi(0)} / (\ell B_U^0)$

Then, $U$ is uniquely defined by $U(\xi) = u(\xi)e^{-z_U \Psi(\xi) - \varepsilon_U \ell \delta \xi}$ for $U = C, P, N$. It
concludes the proof of Lemma 3.1. $\square$

Let us introduce $d(\ell) = \frac{z_C V - \Delta\Psi_1^{pzc} - \Delta\Psi_0^{pzc}}{\ell + \alpha_0 + \alpha_1}$, such that

$$z_C(\Psi(1) - \Psi(0)) = d(\ell)\ell\delta \text{ and } e^{z_C\Psi(\zeta)} = e^{z_C\Psi(0)}e^{d(\ell)\ell\delta\zeta} \quad \forall \zeta \in [0, 1].$$

It permits to rewrite (43a) and (43b) as

(44a) 	$c(\xi) = c(0) - \frac{\hat{J}_C e^{(d(\ell) + \varepsilon_C)\ell\delta\xi} - 1}{\ell\delta} e^{z_C\Psi(0)}$,

(44b) 	$-\frac{\hat{J}_C}{\ell\delta} = \frac{\Gamma_C^1 e^{(d(\ell) + \varepsilon_C)\ell\delta} / B_C^1 - \Gamma_C^0 / B_C^0}{\delta / B_C^0 + \delta e^{(d(\ell) + \varepsilon_C)\ell\delta} / B_C^1 + [e^{(d(\ell) + \varepsilon_C)\ell\delta} - 1] / [d(\ell) + \varepsilon_C]}$.

We are now able to give the existence theorem for (TM2).

**Theorem 3.2.** Let

$$\tilde{K}_1 = \Pi \left( K \gamma_C^1 (\Delta \Psi_1^{pzc}) - \frac{\beta_C^1 (\Delta \Psi_1^{pzc})}{\varepsilon_C} \right) e^{5a_d^0 \tilde{\Psi}},$$

$$\tilde{K}_2 = \Pi K \frac{\beta_C^0 (\tilde{\Psi}) \gamma_C^1 (V - \tilde{\Psi}) - \beta_C^1 (V - \tilde{\Psi}) \gamma_C^0 (\tilde{\Psi})}{\beta_C^0 (\tilde{\Psi}) + \beta_C^1 (V - \tilde{\Psi})} e^{5a_d^0 \tilde{\Psi}},$$

$$\tilde{C} = 1 + \frac{k_d^0 e^{-5a_d^0 \min(\Delta\Psi_0^{pzc}, V - \Delta\Psi_1^{pzc})}}{\Pi K m_C^0 e^{-2b_C^0 \max(\Delta\Psi_0^{pzc}, V - \Delta\Psi_1^{pzc})} + k_C^0 e^{2a_C^0 \min(\Delta\Psi_0^{pzc}, V - \Delta\Psi_1^{pzc})}}.$$

We assume

(45) 	$\min(\tilde{K}_1, \tilde{K}_2) < k_d^0 < \max(\tilde{K}_1, \tilde{K}_2)$.

Then, the problem (TM2) has a solution $(\Psi_s, C_s, \ell_s, \delta_s) \in (C^\infty([0, 1], \mathbb{R}))^2 \times (\mathbb{R}_+^*)^2$
which satisfies $0 \le C_s(\xi) \le \tilde{C}$ for all $\xi \in [0, 1]$.
*Proof.* For a given $\ell$, we will now denote by $(\Psi_\ell, C_\ell, \delta_\ell)$ the solution to (25)-(26)-(27a) given by Lemma 3.1. We introduce the function:

$$F: \quad \mathbb{R}_+^* \to \mathbb{R}$$
$$\ell \mapsto -K \frac{\hat{J}_{C_\ell}}{\ell \delta_\ell}$$

The function $F$ is well defined and continuous. In order to prove the existence of $\ell_s > 0$ such that $F(\ell_s) = 1$, we will apply the intermediate value theorem. Therefore, we have to study the limits of $F$ when $\ell$ tends to 0 and when $\ell$ tends to $+\infty$.

When $\ell \to +\infty$, we have the following limits :

$$\Psi_\ell(0) \to \Delta\Psi_0^{pzc}, \quad \Psi_\ell(1) \to (V - \Delta\Psi_1^{pzc}),$$
$$d(\ell)\ell\delta_\ell \to z_C(V - \Delta\Psi_1^{pzc} - \Delta\Psi_0^{pzc}),$$
$$\delta_\ell \to \frac{k_d^0}{\Pi} e^{-5a_d^0 \Delta\Psi_0^{pzc}}, \quad d(\ell) \to 0.$$

Then, using (44b), we deduce:

$$\lim_{\ell \to +\infty} F(\ell) = K\Pi \frac{\varepsilon_C \gamma_C^1 (\Delta\Psi_1^{pzc})}{\varepsilon_C k_d^0 e^{-5a_d^0 \Delta\Psi_0^{pzc}} + \Pi \beta_C^1 (\Delta\Psi_1^{pzc})}.$$

When $\ell \to 0$, the limits are:

$$\Psi_\ell(0) \to \tilde{\Psi}, \quad \Psi_\ell(1) \to \tilde{\Psi},$$
$$d(\ell)\ell\delta_\ell \to 0, \quad \delta_\ell \to \frac{k_d^0}{\Pi} e^{-5a_d^0 \tilde{\Psi}}.$$

Therefore,

$$\lim_{\ell \to 0} F(\ell) = K \Pi \frac{\gamma_C^1 (V - \tilde{\Psi}) \beta_C^0 (\tilde{\Psi}) - \gamma_C^0 (\tilde{\Psi}) \beta_C^1 (V - \tilde{\Psi})}{k_d^0 e^{-5a_d^0 \tilde{\Psi}} (\beta_C^1 (V - \tilde{\Psi}) + \beta_C^0 (\tilde{\Psi}))}.$$

Then we can see that, if $\min(\tilde{K}_1, \tilde{K}_2) < k_d^0 < \max(\tilde{K}_1, \tilde{K}_2)$,

$$\min \left( \lim_{\ell \to 0} F(\ell), \lim_{\ell \to +\infty} F(\ell) \right) < 1 < \max \left( \lim_{\ell \to 0} F(\ell), \lim_{\ell \to +\infty} F(\ell) \right).$$

It ensures the existence of $\ell_s$ satisfying $F(\ell_s) = 1$. Denoting by $(\Psi_s, C_s, \delta_s)$ the corresponding solution to (25)-(26)-(27a), we obtain that $(\ell_s, \Psi_s, C_s, \delta_s)$ is a solution to (TM2).

Let us now prove that $C_s$ is a positive, bounded function. Writing (43c) for $u = c_s$ while replacing $\hat{J}_C$ by $-\ell_s\delta_s/K$, we obtain that $c_s(0)$ is a sum of positive terms, so that it is positive. Then, writing (44a) for $c = c_s$ and $\xi = 1$, we also get that $c_s(1)$ is positive. Therefore, as $c_s$ is a monotone function thanks to (44a), we conclude that $c_s$ and $C_s$ are positive in $[0, 1]$.
From (42) written for $U = C_s$ and using the known values of $\delta$ and $\hat{J}_C/(\ell\delta)$ at the steady state, we can deduce simple expressions for the boundary values of $C_s$:

$$ \left\{ \begin{array}{l} C_s(0) = \frac{1}{\beta_C^0(\Psi_s(0))} \left( \gamma_C^0(\Psi_s(0)) + \frac{k_d^0 e^{-5a_d^0 \Psi_s(0)}}{\Pi K} \right), \\ C_s(1) = \frac{1}{\beta_C^1(V-\Psi_s(1))} \left( \gamma_C^1(V-\Psi_s(1)) - \frac{k_d^0 e^{-5a_d^0 \Psi_s(0)}}{\Pi K} \right). \end{array} \right. $$

It leads, using (40) and the definition of the functions $\beta_C^0, \beta_C^1, \gamma_C^0, \gamma_C^1$,

$$ C_s(0) \le \frac{\gamma_C^0(\Psi_s(0))}{\beta_C^0(\Psi_s(0))} + \frac{k_d^0 e^{-5a_d^0 \Psi_s(0)}}{\beta_C^0(\Psi_s(0)) \Pi K} \le 1 + \frac{1}{\Pi K} \frac{k_d^0 e^{-5a_d^0 \Psi_s(0)}}{m_C^0 e^{-2b_C^0 \Psi_s(0)} + k_C^0 e^{2a_C^0 \Psi_s(0)}} \le \tilde{C} $$

$$ C_s(1) \le \frac{\gamma_C^1(V - \Psi_s(1))}{\beta_C^1(V - \Psi_s(1))} \le \frac{k_C^1 e^{3a_C^1(V - \Psi_s(1))}}{m_C^1 e^{-3b_C^1(V - \Psi_s(1))} + k_C^1 e^{3a_C^1(V - \Psi_s(1))}} \le 1 \le \tilde{C}. $$

Thanks to the monotonicity of $C_s$, it concludes the proof.

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

**Remark 3.3.** We first remark that the hypothesis (45) reduces to (35) in the case where $\Psi = 0$ (which needs $\Delta\Psi_0^{pzc} = V - \Delta\Psi_1^{pzc} = 0$ and therefore $\tilde{\Psi} = 0$). The hypotheses of Theorem 3.1 and Theorem 3.2 are therefore consistent.

Moreover, Theorem 3.2 is worthy of interest if $\max(\tilde{K}_1, \tilde{K}_2) > 0$. This can be true if either $\beta_C^0(\tilde{\Psi})\gamma_C^1(V - \tilde{\Psi}) \ge \beta_C^1(V - \tilde{\Psi})\gamma_C^0(\tilde{\Psi})$ or $\varepsilon_C K \gamma_C^1(\Delta\Psi_1^{pzc}) \ge \beta_C^1(\Delta\Psi_1^{pzc})$, which is satisfied if

$$ e^{3V} e^{-\tilde{\Psi}} \ge \frac{m_C^0 m_C^1}{k_C^0 k_C^1} \quad \text{or} \quad e^{-3\Delta\Psi_1^{pzc}} \le (\varepsilon_C K - 1) \frac{k_C^1}{m_C^1}. $$

**Remark 3.4.** We do not have any uniqueness result with our method. We can set partial results of uniqueness in the spirit of lemma 3.1. To guarantee the uniqueness of $(\Psi_s, \ell_s, \delta_s, C_s)$ one can impose for instance that $F$ is strictly monotone, which adds new conditions on the parameters.

**Remark 3.5.** For $U = P, N$, there exists a unique solution to (28), given by formula (44a) and (44b). Once again, it is possible to impose conditions such that $\hat{J}_P$ and $\hat{J}_U$ are nonpositive or to ensure that $P_s$ and $N_s$ are non negative.

# 4. NUMERICAL ANALYSIS OF THE SIMPLIFIED MODELS

## 4.1. Presentation of the numerical schemes.
The numerical schemes defined for (TM1) and (TM2) are the adaptation of the scheme (S) introduced in Section 2.2 to these simplified models. They will be respectively denoted (S-TM1)et (S-TM2).
Numerical scheme (S-TM1). The unknowns of the scheme are the discrete densities of oxygen vacancies ($C_i$)$_{0\le i\le I+1}$, the velocity of the interfaces $\delta^h$ and the thickness of the domain $\ell^h$. It writes as follows.

• Scheme for C:

$$ (46a) \qquad \mathcal{G}_{C,i+\frac{1}{2}} - \mathcal{G}_{C,i-\frac{1}{2}} = 0, \quad 0 \le i \le I, $$

$$ (46b) \qquad \mathcal{G}_{C,i+\frac{1}{2}} = \frac{1}{h_{i+\frac{1}{2}}} \left( B\left(h_{i+\frac{1}{2}}(\varepsilon_C \ell^h \delta^h)\right) C_i - B\left(-h_{i+\frac{1}{2}}(\varepsilon_C \ell^h \delta^h)\right) C_{i+1} \right), \quad 1 \le i \le I, $$

$$ (46c) \qquad -\mathcal{G}_{C,\frac{1}{2}} = \ell^h r_C^0(C_0, 0) $$

$$ (46d) \qquad \mathcal{G}_{C,I+\frac{1}{2}} = \ell^h r_C^1(C_{I+1}, 0, V), $$

where B is given by the Bernoulli function (19).

• Scheme for $\delta$ and $\ell$:

$$ (47a) \qquad \delta^h = \frac{k_d^0}{\Pi}, $$

$$ (47b) \qquad \ell^h = -K \frac{\mathcal{G}_{C,I+\frac{1}{2}}}{\delta^h}. $$

Numerical scheme (S-TM2). The unknowns of the scheme are the discrete densities ($C_i$)$_{0\le i\le I+1}$ and the discrete electric potential ($\Psi_i$)$_{0\le i\le I+1}$, the velocity of the interfaces $\delta^h$ and the thickness of the domain $\ell^h$.

• Scheme for $\Psi$:

$$ (48a) \qquad -\frac{\lambda^2}{\ell^h 2} (d\Psi_{i+\frac{1}{2}} - d\Psi_{i-\frac{1}{2}}) = 0, \quad 1 \le i \le I, $$

$$ (48b) \qquad \text{with } d\Psi_{i+\frac{1}{2}} = \frac{\Psi_{i+1} - \Psi_i}{h_{i+\frac{1}{2}}}, \quad 0 \le i \le I, $$

$$ (48c) \qquad \Psi_0 - \frac{\alpha_0}{\ell^h} d\Psi_{\frac{1}{2}} = \Delta\Psi_0^{pzc}, $$

$$ (48d) \qquad \Psi_{I+1} + \frac{\alpha_1}{\ell^h} d\Psi_{I+\frac{1}{2}} = V - \Delta\Psi_1^{pzc}. $$
• Scheme for C:

$$ (49a) \qquad \mathcal{G}_{C,i+\frac{1}{2}} - \mathcal{G}_{C,i-\frac{1}{2}} = 0, \quad 0 \le i \le I, $$

$$ (49b) \qquad \mathcal{G}_{C,i+\frac{1}{2}} = \frac{1}{h_{i+\frac{1}{2}}} \left( B \left( h_{i+\frac{1}{2}} \left( z_C d\Psi_{i+\frac{1}{2}} + \varepsilon_C \ell^h \delta^h \right) \right) C_i - \right. $$

$$ (49c) \qquad \left. B \left( -h_{i+\frac{1}{2}} \left( z_C d\Psi_{i+\frac{1}{2}} + \varepsilon_C \ell^h \delta^h \right) \right) C_{i+1} \right), \quad 1 \le i \le I, $$

$$ (49d) \qquad -\mathcal{G}_{C,\frac{1}{2}} = \ell^h r_C^0(C_0, \Psi_0), $$

$$ (49e) \qquad \mathcal{G}_{C,I+\frac{1}{2}} = \ell^h r_C^1(C_{I+1}, \Psi_{I+1}, V), $$

where B is the Bernoulli function (19).

• Scheme for $\delta$ and $\ell$:

$$ (50a) \qquad \delta^h = \frac{k_d^0}{\Pi} e^{-5a_d^0 \Psi_0}, $$

$$ (50b) \qquad \ell^h = -K \frac{\mathcal{G}_{C,I+\frac{1}{2}}}{\delta^h} $$

In the two next sections, we will show that the schemes (S-TM1) and (S-TM2) are exact. It means that the exact solution to (S-TM1) (respectively (S-TM2)) is a solution to the scheme. It is a well-known property of the Scharfetter-Gummel scheme for a stationary drift-diffusion equation when the drift velocity is the gradient of a potential. The novelty here is that the drift-diffusion equation is coupled with equations on the interface velocity $\delta$ and the size of the layer $\ell$. However, we will prove that the preservation of the steady-state by the Scharfetter-Gummel scheme still holds.

## 4.2. Existence and uniqueness of a solution for (S-TM1).

**Theorem 4.1.** Let assume $[\min(K_1, K_2) < k_d^0 < \max(K_1, K_2)]$, where $K_1$ and $K_2$ are defined by (34). Then, the scheme (S-TM1) has a unique solution $(C^h, \delta^h, \ell^h) \in \mathbb{R}^{I+2} \times (\mathbb{R}_+^*)^2$. Moreover, with $(C_s, \delta_s, \ell_s)$ solution of (TM1), we have $(C^h, \delta^h, \ell^h) = ((C_s(x_i))_{i \in [0, I+1]}, \delta_s, \ell_s)$.

*Proof.* We denote by $(C_s, \delta_s, \ell_s)$ the solution to (TM1) and consider

$$ (C^h, \delta^h, \ell^h) = ((C_s(x_i))_{i \in [0, I+1]}, \delta_s, \ell_s). $$

The formula (30) implies that $(C^h, \delta^h, \ell^h)$ fulfills (46) where, for all $i \in [0, I]$, $\mathcal{G}_{C^h, i+\frac{1}{2}} = \hat{J}_{C_s}$. Then (47) is exactly (23).
In order to prove the uniqueness, let us consider $(C^h, \delta^h, \ell^h)$ solution of (S-TM1). Equations (47) imply

$$ \delta^h = \frac{k_d^0}{\Pi} = \delta_s, \quad \frac{\mathcal{G}_{C^h}}{\ell^h} = -\frac{k_d^0}{K\Pi} = \frac{\hat{J}_{C_s}}{\ell_s}, $$

$$ C_0^h = \frac{1}{\beta_C^0(0)} \left( \gamma_C^0(0) + \frac{k_d^0}{\Pi K} \right) = C_s(0), $$

$$ C_1^h = \frac{1}{\beta_C^1(V)} \left( \gamma_C^1(V) - \frac{k_d^0}{\Pi K} \right) = C_s(1). $$

Moreover by definition of $C^h$ and the Bernoulli function we have for any $i \in [0, I]$

$$ (51) \qquad C_{i+1}^h = \left( C_i^h + \frac{\mathcal{G}_{C^h}}{\varepsilon_C \ell^h \delta_s} \right) e^{-\varepsilon_C \ell^h \delta_s (x_{i+1} - x_i)} - \frac{\mathcal{G}_{C^h}}{\varepsilon_C \ell^h \delta_s}. $$

Starting at $i = I$ and going down to $i = 0$ leads to

$$ \begin{aligned} C_{I+1}^h &= -\frac{\mathcal{G}_{C^h}}{\varepsilon_C \ell^h \delta_s} + \left( \frac{\mathcal{G}_{C^h}}{\varepsilon_C \ell^h \delta_s} + C_I^h \right) e^{-\frac{\ell^h \delta_s}{D} h_{I+1/2}} \\ &= -\frac{\mathcal{G}_{C^h}}{\varepsilon_C \ell^h \delta_s} + \left( \frac{\mathcal{G}_{C^h}}{\varepsilon_C \ell^h \delta_s} + C_{I-1}^h \right) e^{-\frac{\ell^h \delta_s}{D} (h_{I+1/2} + h_{I-1/2})} \\ &= -\frac{\mathcal{G}_{C^h}}{\varepsilon_C \ell^h \delta_s} + \left( \frac{\mathcal{G}_{C^h}}{\varepsilon_C \ell^h \delta_s} + C_0^h \right) e^{-\frac{\ell^h \delta_s}{D}}. \end{aligned} $$

Since $C_0^h = C_s(0)$, $C_{I+1}^h = C_s(1)$ and $\frac{\mathcal{G}_{C^h}}{\ell^h \delta_s} = \frac{\hat{J}_{C_s}}{\ell_s \delta_s}$ we deduce, thanks to (32), that

$$ \ell^h = \frac{\Pi}{\varepsilon_C k_d^0} \ln \left( \frac{\varepsilon_C C_0^h - 1/K}{\varepsilon_C C_1^h - 1/K} \right) = \ell_s. $$

Therefore, $\mathcal{G}_{C^h} = \hat{J}_{C_s}$ and, using (51), for any $i \in [0, I+1]$, $C_i^h = C_s(x_i)$. It concludes the proof. $\square$

### 4.3. Existence of a solution for (S-TM2).

**Theorem 4.2.** We assume

$$ \min (\tilde{K}_1, \tilde{K}_2) < k_d^0 < \max (\tilde{K}_1, \tilde{K}_2), $$

where $\tilde{K}_1, \tilde{K}_2$ are defined in Theorem 3.2. Then the scheme (S-TM2) has a solution $(\Psi^h, C^h, \delta^h, \ell^h) \in (\mathbb{R}^{I+2})^2 \times (\mathbb{R}_+^*)^2$. Moreover there exists $(\Psi_s, C_s, \delta_s, \ell_s)$ solution to (TM2) such that:

$$ (52) \qquad \forall i \in [0, I+1], \quad (\Psi_i^h, C_i^h, \delta^h, \ell^h) = (\Psi_s(ih), C_s(ih), \delta_s, \ell_s). $$
*Proof. Existence.* Let $(\Psi_s, C_s, \delta_s, \ell_s)$ be a solution of (TM2). Since (48) is exact on the linear solutions of (25), we define $\Psi_i^h = \Psi(ih)$, $i \in [0, I+1]$ solution of (48) with $\ell^h = \ell_s$. We then fix

$$ \delta^h = \delta_s = \frac{k_d^0}{\Pi} e^{-5a_d^0 \Psi_0^h}. $$

We observe that $(z_C d\Psi_{i+\frac{1}{2}} + \varepsilon_C \ell^h \delta^h)$ is constant, therefore the Scharfetter-Gummel scheme (49) is exact on (26). Define $(C_i^h) = (C(ih))$, $i \in [0, I+1]$, it is solution of (49). Finally (27b) implies (50b). It concludes the proof of the existence of a solution.

*Any solution is exact.* Let $(\Psi^h, C^h, \delta^h, \ell^h)$ be a solution of (S-TM2). We apply lemma 3.1 with $\ell = \ell^h$ it gives $(\Psi_{\ell^h}, C_{\ell^h}, \delta_{\ell^h}) \in (C^\infty[0, 1])^2 \times \mathbb{R}_+^*$, solution of (25)-(26)-(27a). Since $\Psi_{\ell^h}$ is linear and the scheme given by (48) is exact on linear function we have for all $i \in [0, I+1]$ $\Psi_i^h = \Psi_{\ell^h}(ih)$. Therefore (50a) implies $\delta^h = \delta_{\ell^h}$. The exactness of the Scharfetter-Gummel scheme gives $C_i^h = C_{\ell^h}(ih)$, for all $i \in [0, I+1]$. Then equations (49e) and (26c) implies $G_{C^h, I+\frac{1}{2}} = \hat{J}_{C_{\ell^h}}$. Combined to (50b) we deduce that

$$ \ell^h = -K \frac{G_{C^h, I+\frac{1}{2}}}{\delta^h} = -K \frac{\hat{J}_{C_{\ell^h}}}{\delta_{\ell^h}}. $$

It shows that $(\Psi_{\ell^h}, C_{\ell^h}, \delta_{\ell^h}, \ell^h)$ is a solution to (TM2) and concludes the proof.

□

## APPENDIX A. PARAMETER

All the numerical simulations have been done with the following set of parameters. The scaling process and the physical values used to obtain the scaled values in A.2 are detailed in section 5 of [6].

<table><caption>TABLE A.1. Physical parameters of the test case</caption><thead><tr><th>D<sub>1</sub> (m<sup>2</sup> &middot; s<sup>-1</sup>)</th><th>D<sub>2</sub> (m<sup>2</sup> &middot; s<sup>-1</sup>)</th><th>D<sub>3</sub> (m<sup>2</sup> &middot; s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>10<sup>-20</sup></td><td>10<sup>-6</sup></td><td>10<sup>-20</sup></td></tr><tr><td>m</td><td>&Omega;<sub>ox</sub> (m<sup>3</sup> &middot; mol)</td><td>&Omega;<sub>Fe</sub> (m<sup>3</sup> &middot; mol)</td></tr><tr><td>3</td><td>4.474 &middot; 10<sup>-5</sup></td><td>7.105 &middot; 10<sup>-6</sup></td></tr><tr><td>(a<sub>u</sub><sup>0</sup>, b<sub>u</sub><sup>0</sup>) (u = P, N, C, r)</td><td>(a<sub>u</sub><sup>1</sup>, b<sub>u</sub><sup>1</sup>) (u = P, C)</td><td>a<sub>d</sub><sup>0</sup></td></tr><tr><td>(0.5, 0.5)</td><td>(0.5, 0.5)</td><td>0</td></tr></tbody></table>
<table><tr><td>ΔΨ<sup>pzc</sup><sub>0</sub></td><td>ΔΨ<sup>pzc</sup><sub>1</sub></td><td>k<sup>0</sup><sub>d,ref</sub></td></tr><tr><td>-8.22136 · 10<sup>-12</sup></td><td>0</td><td>64.267 · 10<sup>-pH</sup></td></tr><tr><td>P<sup>m</sup></td><td>N<sub>metal</sub></td><td>ρ<sub>hl</sub></td></tr><tr><td>2</td><td>2.48502 · 10<sup>-1</sup></td><td>-5</td></tr><tr><td>λ<sup>2</sup></td><td>α<sub>0</sub></td><td>α<sub>1</sub></td></tr><tr><td>1.05432 · 10<sup>-3</sup></td><td>0.177083</td><td>8.854 · 10<sup>-2</sup></td></tr></table>

<table><caption>TABLE A.2. Scaled parameters of the test case</caption><tr><td>m<sup>0</sup><sub>P</sub></td><td>k<sup>0</sup><sub>P</sub></td><td>m<sup>1</sup><sub>P</sub></td><td>k<sup>1</sup><sub>P</sub></td></tr><tr><td>0</td><td>10<sup>8</sup></td><td>10<sup>8</sup></td><td>10<sup>14</sup></td></tr><tr><td>m<sup>0</sup><sub>N</sub></td><td>k<sup>0</sup><sub>N</sub></td><td>m<sup>1</sup><sub>N</sub></td><td>k<sup>1</sup><sub>N</sub></td></tr><tr><td>0</td><td>2.10740 · 10<sup>-19</sup></td><td>2.68111 · 10<sup>1</sup></td><td>2.68111 · 10<sup>1</sup></td></tr><tr><td>m<sup>0</sup><sub>C</sub></td><td>k<sup>0</sup><sub>C</sub></td><td>m<sup>1</sup><sub>C</sub></td><td>k<sup>1</sup><sub>C</sub></td></tr><tr><td>1.78113 · 10<sup>-13</sup></td><td>4.474 · 10<sup>47</sup></td><td>0</td><td>7.85325</td></tr></table>

REFERENCES

[1] T. Aiki and A. Muntean. Existence and uniqueness of solutions to a mathematical model predicting service life of concrete structures. *Adv. Math. Sci. Appl.*, 19(1):109–129, 2009.

[2] T. Aiki and A. Muntean. On uniqueness of a weak solution of one-dimensional concrete carbonation problem. *Discrete Contin. Dyn. Syst.*, 29(4):1345–1365, 2011.

[3] T. Aiki and A. Muntean. A free-boundary problem for concrete carbonation: front nucleation and rigorous justification of the $\sqrt{t}$-law of propagation. *Interfaces Free Bound.*, 15(2):167–180, 2013.

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

[5] C. Bataillon, F. Bouchon, C. Chainais-Hillairet, J. Fuhrmann, E. Hoarau, and R. Touzani. Numerical methods for simulation of a corrosion model with moving numerical methods for simulation of a corrosion model with moving oxide layer. *Journal of Computational Physics*, 231(18):6213–6231, 2012.

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

[7] C. Chainais-Hillairet, P.-L. Colin, and I. Lacroix-Violet. Convergence of a Finite Volume Scheme for a Corrosion Model. *Sept. 2014*.

[8] C. Chainais-Hillairet and I. Lacroix-Violet. The existence of solutions to a corrosion model. *Appl. Math. Lett.*, 25(11):1784–1789, 2012.

[9] C. Chainais-Hillairet and I. Lacroix-Violet. On the existence of solutions for a driftdiffusion system arising in corrosion modelling. *DCDS-B*, 2014.
[10] 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.

[11] H. Gajewski and K. Gärtner. On the discretization of van Roosbroeck's equations with magnetic field. *Z. Angew. Math. Mech.*, 76(5):247–264, 1996.

[12] A. M. Il'in. A difference scheme for a differential equation with a small parameter multiplying the highest derivative. *Mat. Zametki*, 6:237–248, 1969.

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

[14] P. A. Markowich, C. A. Ringhofer, and C. Schmeiser. *Semiconductor equations*. Springer-Verlag, Vienna, 1990.

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

UNIVERSITÉ DE LILLE, LABORATOIRE PAUL PAINLEVÉ, CITÉ SCIENTIFIQUE, 59655 VILLENEUVE D'ASCQ CEDEX, FRANCE

E-mail address: Claire.Chainais@math.univ-lille1.fr

UNIVERSITÉ LIBRE DE BRUXELLES (ULB), DÉPARTEMENT DE MATHÉMATIQUE, BOULEVARD DU TRIOMPHE 1050 BRUSSELS, BELGIUM AND PROJECT-TEAM MEPHYSTO, INRIA LILLE - NORD EUROPE, 59655 VILLENEUVE D'ASCQ, FRANCE

E-mail address: Thomas.Gallouet@ulb.ac.be