# Anodic dissolution model with diffusion-migration transport for simulating localized corrosion

Meriem Bouguezzi, Jean-François Scheid, Danielle Hilhorst, Hiroshi Matano,
Christian Bataillon, Florence Lequien, Fabien Rouillard

► To cite this version:

Meriem Bouguezzi, Jean-Francois Scheid, Danielle Hilhorst, Hiroshi Matano, Christian Bataillon, et al. Anodic dissolution model with diffusion-migration transport for simulating localized corrosion. Electrochimica Acta, 2024, 477, pp.143806. 10.1016/j.electacta.2024.143806. hal-04160934v2

HAL Id: hal-04160934

https://inria.hal.science/hal-04160934v2

Submitted on 9 Jan 2024

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
# Anodic dissolution model with diffusion-migration transport for simulating localized corrosion

Meriem Bouguezzi, Jean-François Scheid, Danielle Hilhorst, Hiroshi Matano,
Christian Bataillon, Florence Lequien, Fabien Rouillard

► To cite this version:

Meriem Bouguezzi, Jean-Francois Scheid, Danielle Hilhorst, Hiroshi Matano, Christian Bataillon, et al. Anodic dissolution model with diffusion-migration transport for simulating localized corrosion. Electrochimica Acta, 2024, 477, pp.143806. 10.1016/j.electacta.2024.143806. hal-04160934v2

HAL Id: hal-04160934

https://inria.hal.science/hal-04160934v2

Submitted on 9 Jan 2024

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
# Anodic dissolution model with diffusion-migration transport for simulating localized corrosion

Meriem Bouguezzi<sup>a,e</sup>, Jean-François Scheid<sup>b,*</sup>, Danielle Hilhorst<sup>c</sup>, Hiroshi Matano<sup>d</sup>, Christian Bataillon<sup>e</sup>, Florence Lequien<sup>e</sup>, Fabien Rouillard<sup>e</sup>

<sup>a</sup>Laboratoire de Mathématiques d'Orsay, Université Paris-Saclay, 91405, Orsay Cedex, France

<sup>b</sup>Université de Lorraine, CNRS, Inria, IECL, Nancy, France

<sup>c</sup>CNRS and Laboratoire de Mathématiques d'Orsay, Université Paris-Saclay, 91405, Orsay Cedex, France

<sup>d</sup>Meiji Institute for Advanced Study of Mathematical Sciences, Meiji University, 4-21-1 Nakano, 164-8525, Tokyo, Japan

<sup>e</sup>Université Paris-Saclay, CEA, Service de recherche en Corrosion et Comportement des Matériaux, 91191, Gif-sur-Yvette, France

## Abstract

A one-dimensional model for simulating localized corrosion is developed, incorporating a time-dependent diffusion-migration transport and a moving interface. The model is applied to iron with the Butler-Volmer formula as the dissolution law, and both crevice and pit configurations are simulated. The mathematical procedure for solving this strongly coupled differential equations system and the numerical development for simulations are detailed. A finite-difference ALE scheme is used for the numerical computation of the solutions of this free boundary problem. The results show that the dissolution rate increases with chloride concentration and metal potential, and migration plays, most of the time, a significant role in species transport for both crevice and pit configurations. The evolution of the repassivation potential with pit depth is computed and is in good agreement with experimental results. The time-dependent model under consideration provides more accurate quantitative results for concentration profile and crevice depth than stationary models.

**Keywords:** Stefan problem, free boundary problem, localized corrosion, diffusion-migration transport, finite-difference ALE scheme.

## 1. Introduction

One of the major causes of structural material failure in aggressive aqueous solutions containing chlorides is localized corrosion, such as crevice or pitting corrosion. The crevice corrosion occurs within narrow gaps or crevices of metallic components that already exist before the corrosion starts. On the other hand, the pitting corrosion is initiated typically by an imperfection on an initially flat metallic surface. Therefore, in this type of corrosion, the size of the pit is initially very small, though it may eventually grow much larger under certain circumstances. In either case, such corrosion processes can result in localized failure in a very short period of time. Many localized corrosion models, mainly for crevice corrosion, have been proposed in literature to study and predict the time evolution of localized corrosion propagation. In the literature, there are two typical approaches for building mathematical models for corrosion: empirical modeling and mechanistic based modeling, the latter being based on the fundamental laws of mass transport, chemistry

<footnote>*Corresponding author</footnote>
# Anodic dissolution model with diffusion-migration transport for simulating localized corrosion

Meriem Bouguezzi<sup>a,e</sup>, Jean-François Scheid<sup>b,*</sup>, Danielle Hilhorst<sup>c</sup>, Hiroshi Matano<sup>d</sup>, Christian Bataillon<sup>e</sup>, Florence Lequien<sup>e</sup>, Fabien Rouillard<sup>e</sup>

<sup>a</sup>Laboratoire de Mathématiques d'Orsay, Université Paris-Saclay, 91405, Orsay Cedex, France

<sup>b</sup>Université de Lorraine, CNRS, Inria, IECL, Nancy, France

<sup>c</sup>CNRS and Laboratoire de Mathématiques d'Orsay, Université Paris-Saclay, 91405, Orsay Cedex, France

<sup>d</sup>Meiji Institute for Advanced Study of Mathematical Sciences, Meiji University, 4-21-1 Nakano, 164-8525, Tokyo, Japan

<sup>e</sup>Université Paris-Saclay, CEA, Service de recherche en Corrosion et Comportement des Matériaux, 91191, Gif-sur-Yvette, France

## Abstract

A one-dimensional model for simulating localized corrosion is developed, incorporating a time-dependent diffusion-migration transport and a moving interface. The model is applied to iron with the Butler-Volmer formula as the dissolution law, and both crevice and pit configurations are simulated. The mathematical procedure for solving this strongly coupled differential equations system and the numerical development for simulations are detailed. A finite-difference ALE scheme is used for the numerical computation of the solutions of this free boundary problem. The results show that the dissolution rate increases with chloride concentration and metal potential, and migration plays, most of the time, a significant role in species transport for both crevice and pit configurations. The evolution of the repassivation potential with pit depth is computed and is in good agreement with experimental results. The time-dependent model under consideration provides more accurate quantitative results for concentration profile and crevice depth than stationary models.

**Keywords:** Stefan problem, free boundary problem, localized corrosion, diffusion-migration transport, finite-difference ALE scheme.

## 1. Introduction

One of the major causes of structural material failure in aggressive aqueous solutions containing chlorides is localized corrosion, such as crevice or pitting corrosion. The crevice corrosion occurs within narrow gaps or crevices of metallic components that already exist before the corrosion starts. On the other hand, the pitting corrosion is initiated typically by an imperfection on an initially flat metallic surface. Therefore, in this type of corrosion, the size of the pit is initially very small, though it may eventually grow much larger under certain circumstances. In either case, such corrosion processes can result in localized failure in a very short period of time. Many localized corrosion models, mainly for crevice corrosion, have been proposed in literature to study and predict the time evolution of localized corrosion propagation. In the literature, there are two typical approaches for building mathematical models for corrosion: empirical modeling and mechanistic based modeling, the latter being based on the fundamental laws of mass transport, chemistry

<footnote>*Corresponding author</footnote>
and electrochemistry processes. Both approaches have their advantages and disadvantages, but the mechanistic based models, which we discuss in the present article, are considered to be more flexible [1]. Localized corrosion involves the dissolution of metal, mass transport by diffusion and migration, and reactions in solution. From a mathematical point of view, this problem can be identified as a Stefan problem involving a convection-reaction-diffusion system. The mechanistic based models proposed in literature differ from each other in the considered mass transport equation (diffusion with or without migration, with or without reactions between species), the coupling or not of the mass transport equation with a moving boundary, and the equation used for the electric potential field (Poisson equation or electroneutrality condition). These differences have a strong impact on the mathematical expression of the resulting coupled differential equations. Moreover, different boundary conditions and initial conditions have been used, mainly in 1 or 2 space dimensions. Finally, several computational methodologies have been proposed to solve the coupled differential equations and to follow the movement of the free boundary, when considered. The displacement of the moving interface may be tackled by ALE methods with moving mesh techniques (see [2] for a moving vertical crevice wall model and also [3] for a crevice geometry with axial symmetry). A phase-field approach may also be considered as in [4], leading to a different Partial Differential Equation System from sharp interface models. In some studies, the moving interface is maintained fixed ([5],[6],[7]) or reduced models are considered [8]. We also mention [9] (see also [10]) where a rigorous derivation of a simple model of pit propagation from the mass conservation and Fick's law is presented with a 2D finite volume method of resolution.

In many studies, finite element commercial softwares such as COMSOL are used. Even if all the results obtained by these models are quite convincing, they are, for most of them, difficult to reproduce, without using these softwares, since the mathematical solution for the coupled differential equations and the computational methodology with the numerical scheme are not completely described. A detailed description of the mathematical equations used for simulating crevice corrosion and of the assumptions often made in many papers to simplify their resolution have been recently given by Moraes et al. [11]. Nevertheless, this paper gives few details on the mathematical and numerical procedures used to solve the coupled differential equation system.

The main objective of this paper is to write and justify as clearly as possible the mathematical description of a localized corrosion problem and to explain in detail, step-by-step, how the coupled differential equation system is solved by a well suited numerical scheme. Our objective is to provide corrosion specialists with a well-detailed and justified mathematical description, followed by a robust computational methodology, so that corrosion specialists can develop their own numerical model without using any commercial software.

The localized corrosion problem described in this paper is the dissolution of pure iron in aqueous solutions containing NaCl, which is the simplest physical system for this corrosion type. A time-dependent localized corrosion model, taking into account mass transport by diffusion and migration, and dissolution of iron according to the Butler-Volmer law, is proposed, and its mathematical structure and the method for solving the problem numerically are described in detail (see appendices). The consideration of the Butler-Volmer law allows to simulate dissolution at any metal potentials and not only far from equilibrium as it is the case when Tafel's law is used.

First, the physical model governing pit propagation is described in detail. Then, the coupled differential equations derived from such a model are given, and a mathematical method for solving the equations is proposed. Finally, a numerical scheme is described as accurately as possible to help the understanding of non-specialists.

The numerical scheme has been implemented in an open-source software (Python). Several numerical simulations are presented as a function of localized configuration, crevice, or pit. First,
and electrochemistry processes. Both approaches have their advantages and disadvantages, but the mechanistic based models, which we discuss in the present article, are considered to be more flexible [1]. Localized corrosion involves the dissolution of metal, mass transport by diffusion and migration, and reactions in solution. From a mathematical point of view, this problem can be identified as a Stefan problem involving a convection-reaction-diffusion system. The mechanistic based models proposed in literature differ from each other in the considered mass transport equation (diffusion with or without migration, with or without reactions between species), the coupling or not of the mass transport equation with a moving boundary, and the equation used for the electric potential field (Poisson equation or electroneutrality condition). These differences have a strong impact on the mathematical expression of the resulting coupled differential equations. Moreover, different boundary conditions and initial conditions have been used, mainly in 1 or 2 space dimensions. Finally, several computational methodologies have been proposed to solve the coupled differential equations and to follow the movement of the free boundary, when considered. The displacement of the moving interface may be tackled by ALE methods with moving mesh techniques (see [2] for a moving vertical crevice wall model and also [3] for a crevice geometry with axial symmetry). A phase-field approach may also be considered as in [4], leading to a different Partial Differential Equation System from sharp interface models. In some studies, the moving interface is maintained fixed ([5],[6],[7]) or reduced models are considered [8]. We also mention [9] (see also [10]) where a rigorous derivation of a simple model of pit propagation from the mass conservation and Fick's law is presented with a 2D finite volume method of resolution.

In many studies, finite element commercial softwares such as COMSOL are used. Even if all the results obtained by these models are quite convincing, they are, for most of them, difficult to reproduce, without using these softwares, since the mathematical solution for the coupled differential equations and the computational methodology with the numerical scheme are not completely described. A detailed description of the mathematical equations used for simulating crevice corrosion and of the assumptions often made in many papers to simplify their resolution have been recently given by Moraes et al. [11]. Nevertheless, this paper gives few details on the mathematical and numerical procedures used to solve the coupled differential equation system.

The main objective of this paper is to write and justify as clearly as possible the mathematical description of a localized corrosion problem and to explain in detail, step-by-step, how the coupled differential equation system is solved by a well suited numerical scheme. Our objective is to provide corrosion specialists with a well-detailed and justified mathematical description, followed by a robust computational methodology, so that corrosion specialists can develop their own numerical model without using any commercial software.

The localized corrosion problem described in this paper is the dissolution of pure iron in aqueous solutions containing NaCl, which is the simplest physical system for this corrosion type. A time-dependent localized corrosion model, taking into account mass transport by diffusion and migration, and dissolution of iron according to the Butler-Volmer law, is proposed, and its mathematical structure and the method for solving the problem numerically are described in detail (see appendices). The consideration of the Butler-Volmer law allows to simulate dissolution at any metal potentials and not only far from equilibrium as it is the case when Tafel's law is used.

First, the physical model governing pit propagation is described in detail. Then, the coupled differential equations derived from such a model are given, and a mathematical method for solving the equations is proposed. Finally, a numerical scheme is described as accurately as possible to help the understanding of non-specialists.

The numerical scheme has been implemented in an open-source software (Python). Several numerical simulations are presented as a function of localized configuration, crevice, or pit. First,
the simulated results for crevice configuration are compared to results from the Sharland model [6, 7], and the influence of physical parameters such as the applied electrical potential on iron and the chloride concentration in the solution on the corrosion rate is analyzed. Then, the proportion of diffusion and migration in the transport of the species in crevice and pit configurations is calculated as a function of chloride concentration and metal potential. Finally, the evolution of the repassivation potential with pit depth is calculated from the model. The results are in good agreement with the experimental observations.

## 2. Diffusion-migration model with moving interface

As we mentioned earlier, there are typically two types of localized corrosions: crevice corrosion and pitting corrosion. Crevice corrosion occurs in a narrow gap or crevice in a metallic body that traps corrosive agents. Therefore, the initial size of the crevice is relatively large. In contrast, pitting corrosion usually occurs on a flat surface of the metallic body when there is some sort of imperfection on the surface. Therefore, the initial size of the pit for the pitting corrosion is tiny in general. When pitting corrosion starts, there often appears a porous metallic cap, or a lacy cover, that covers the mouth of the pit, and a salt film with high concentration of chloride ions at the bottom of the pit. The presence of a metallic cover accelerates corrosion inside the pit. When the metallic cover disappears, the speed of the corrosion slows down sharply, but by that time the pit may have grown to a good size. In this paper, in addition to crevice corrosion, we deal with the pitting corrosion after the metallic cover has disappeared. We are interested in investigating the conditions that allow the pit to grow without the presence of the metallic cover. The effect of a salt film at the bottom will be taken into account. These two types of localized corrosions can be treated by the same mathematical model, though one has to choose different parameter values for each type of the corrosion; see subsections 4.3 and 4.4.

### 2.1. *Derivation of the equations*

We consider a simple situation where corrosion only occurs at the bottom of the pit/crevice and the pit/crevice has a thin cylindrical shape (see Figure 2.1). We assume that the pit/crevice is so thin that the aqueous solution is well mixed by diffusion in the direction orthogonal to the axis of the cylinder. Note that a one-dimensional pit should be represented as an interval with a moving boundary $[0, x_d(t)]$ but in order to show the physical context, we represent it as a rectangle in Figure 2.1. In such a situation, it is reasonable to assume that the concentration of the ions and the potential $\phi$ depend on a one space variable $x$, so that the system reduces to a one-dimensional problem.
the simulated results for crevice configuration are compared to results from the Sharland model [6, 7], and the influence of physical parameters such as the applied electrical potential on iron and the chloride concentration in the solution on the corrosion rate is analyzed. Then, the proportion of diffusion and migration in the transport of the species in crevice and pit configurations is calculated as a function of chloride concentration and metal potential. Finally, the evolution of the repassivation potential with pit depth is calculated from the model. The results are in good agreement with the experimental observations.

## 2. Diffusion-migration model with moving interface

As we mentioned earlier, there are typically two types of localized corrosions: crevice corrosion and pitting corrosion. Crevice corrosion occurs in a narrow gap or crevice in a metallic body that traps corrosive agents. Therefore, the initial size of the crevice is relatively large. In contrast, pitting corrosion usually occurs on a flat surface of the metallic body when there is some sort of imperfection on the surface. Therefore, the initial size of the pit for the pitting corrosion is tiny in general. When pitting corrosion starts, there often appears a porous metallic cap, or a lacy cover, that covers the mouth of the pit, and a salt film with high concentration of chloride ions at the bottom of the pit. The presence of a metallic cover accelerates corrosion inside the pit. When the metallic cover disappears, the speed of the corrosion slows down sharply, but by that time the pit may have grown to a good size. In this paper, in addition to crevice corrosion, we deal with the pitting corrosion after the metallic cover has disappeared. We are interested in investigating the conditions that allow the pit to grow without the presence of the metallic cover. The effect of a salt film at the bottom will be taken into account. These two types of localized corrosions can be treated by the same mathematical model, though one has to choose different parameter values for each type of the corrosion; see subsections 4.3 and 4.4.

### 2.1. *Derivation of the equations*

We consider a simple situation where corrosion only occurs at the bottom of the pit/crevice and the pit/crevice has a thin cylindrical shape (see Figure 2.1). We assume that the pit/crevice is so thin that the aqueous solution is well mixed by diffusion in the direction orthogonal to the axis of the cylinder. Note that a one-dimensional pit should be represented as an interval with a moving boundary $[0, x_d(t)]$ but in order to show the physical context, we represent it as a rectangle in Figure 2.1. In such a situation, it is reasonable to assume that the concentration of the ions and the potential $\phi$ depend on a one space variable $x$, so that the system reduces to a one-dimensional problem.
<figure><img src="image_2.png" /><figcaption>Figure 2.1: A one dimensional corrosion pit.</figcaption></figure>

It is composed of the following domains :

1. Solid metal with temporally and spatially constant metal atom concentration.

2. Pit solution with temporally and spatially varying ion concentrations.

3. Bulk solution outside of the pit/crevice, characterized by the vanishing iron ion concentration $C_{Fe^{2+}} = 0$.

The depth of the pit/crevice is given by $x_d(t)$ which indicates the position of the moving boundary at time $t > 0$ and the entrance of the pit/crevice corresponds to the position $x = 0$.

We suppose that the pit/crevice solution contains $M$ species of ions, whose concentrations are denoted by $C_j = C_j(x, t)$ ($j = 1, \dots, M$) at the position $x$ and time $t$. The concentrations $C_j$ and the electrostatic potential $\phi = \phi(x, t)$ satisfy the following Poisson-Nernst-Planck system [1] :

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right), \quad t > 0, 0 < x < x_d(t), \quad j = 1, \dots, M, \quad (2.1a) $$

$$ -\frac{\partial}{\partial x} \left( \epsilon \frac{\partial \phi}{\partial x} \right) = \sum_{j=1}^{M} z_j F C_j, \quad t > 0, 0 < x < x_d(t), \quad (2.1b) $$

where $D_j, z_j$ denote respectively, the diffusion coefficient and the valence of the ion species $j$, $F$ the Faraday constant, $R$ the ideal gas constant, $T$ the temperature and $\epsilon$ the dielectric constant [12, p.96].

On one hand, the time dependent form of the Nernst-Planck equation (2.1a) is a conservation of mass equation used to describe the motion of charged chemical species in the aqueous pit/crevice solution $[0, x_d(t))$. It extends the diffusion Fick's law to the case where the diffusing particles are also moved by electrostatic forces. These electrostatic forces are caused by the electric field $E = -\frac{\partial \phi}{\partial x}$. On the other hand, the Poisson equation (2.1b) describes the electrical potential field $E$ which results from the charge distribution in the pit solution.

As the dielectric constant $\epsilon$ is extremely small in general, the left-hand side of (2.1b) is negligible under the typical scaling of the phenomena that we discuss. Therefore the Poisson equation (2.1b) can be replaced by the local electroneutrality condition which is given by

$$ \sum_{j=1}^{M} z_j C_j = 0. \quad (2.2) $$
<figure><img src="image_2.png" /><figcaption>Figure 2.1: A one dimensional corrosion pit.</figcaption></figure>

It is composed of the following domains :

1. Solid metal with temporally and spatially constant metal atom concentration.

2. Pit solution with temporally and spatially varying ion concentrations.

3. Bulk solution outside of the pit/crevice, characterized by the vanishing iron ion concentration $C_{Fe^{2+}} = 0$.

The depth of the pit/crevice is given by $x_d(t)$ which indicates the position of the moving boundary at time $t > 0$ and the entrance of the pit/crevice corresponds to the position $x = 0$.

We suppose that the pit/crevice solution contains $M$ species of ions, whose concentrations are denoted by $C_j = C_j(x, t)$ ($j = 1, \dots, M$) at the position $x$ and time $t$. The concentrations $C_j$ and the electrostatic potential $\phi = \phi(x, t)$ satisfy the following Poisson-Nernst-Planck system [1] :

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right), \quad t > 0, 0 < x < x_d(t), \quad j = 1, \dots, M, \quad (2.1a) $$

$$ -\frac{\partial}{\partial x} \left( \epsilon \frac{\partial \phi}{\partial x} \right) = \sum_{j=1}^{M} z_j F C_j, \quad t > 0, 0 < x < x_d(t), \quad (2.1b) $$

where $D_j, z_j$ denote respectively, the diffusion coefficient and the valence of the ion species $j$, $F$ the Faraday constant, $R$ the ideal gas constant, $T$ the temperature and $\epsilon$ the dielectric constant [12, p.96].

On one hand, the time dependent form of the Nernst-Planck equation (2.1a) is a conservation of mass equation used to describe the motion of charged chemical species in the aqueous pit/crevice solution $[0, x_d(t))$. It extends the diffusion Fick's law to the case where the diffusing particles are also moved by electrostatic forces. These electrostatic forces are caused by the electric field $E = -\frac{\partial \phi}{\partial x}$. On the other hand, the Poisson equation (2.1b) describes the electrical potential field $E$ which results from the charge distribution in the pit solution.

As the dielectric constant $\epsilon$ is extremely small in general, the left-hand side of (2.1b) is negligible under the typical scaling of the phenomena that we discuss. Therefore the Poisson equation (2.1b) can be replaced by the local electroneutrality condition which is given by

$$ \sum_{j=1}^{M} z_j C_j = 0. \quad (2.2) $$
In what follows we shall always assume that the solution is completely electrically neutral in the entire pit/crevice solution.

Next we multiply the equations in (2.1a) by $z_j F$ and sum them up for $j = 1, \dots, M$. With equation (2.2), this leads to

$$z_1 F \frac{\partial C_j}{\partial t} + \dots + z_M F \frac{\partial C_M}{\partial t} = F \frac{\partial}{\partial t} (z_1 C_1 + \dots + z_M C_M) = 0.$$

Consequently, we obtain from equation (2.1a),

$$\sum_{j=1}^{M} \frac{\partial}{\partial x} \left( D_j \left( z_j F \frac{\partial C_j}{\partial x} + \frac{z_j^2 F^2}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right) = 0.$$

This equation can be rewritten as

$$\frac{\partial}{\partial x} \left( \sigma \frac{\partial \phi}{\partial x} \right) + \sum_{j=1}^{M} \frac{\partial}{\partial x} \left( z_j F D_j \frac{\partial C_j}{\partial x} \right) = 0, \quad t > 0, 0 < x < x_d(t) \quad \text{where} \quad \sigma := \sum_{j=1}^{M} \frac{z_j^2 F^2}{RT} D_j C_j, \quad (2.3)$$

and $\sigma$ is the electrical conductivity of the ionic solution.

Thus, in order to keep the electroneutrality of the solution, $\phi$ has to satisfy (2.3). In the special case where $C_j$ are spatially uniform (i.e. $\frac{\partial C_j}{\partial x} = 0, j = 1, \dots, M$), equation (2.3) reduces to Laplace's equation $\frac{\partial^2 \phi}{\partial x^2} = 0$, hence $\phi$ is harmonic. However, $\phi$ is not harmonic except in such special cases, and thus in general, the Laplace equation cannot be used (see [13] for a general discussion on this topic).

## 2.2. Initial and boundary conditions

Henceforth, the equation system (2.1)–(2.2) is considered. Needless to say, (2.1)–(2.2) cannot be solved without specifying the boundary conditions as well as the initial condition at $t = 0$. Since the domain $[0, x_d(t))$, which represents the pit/crevice, evolves along with $\{C_j\}$ and $\phi$, this is a free boundary problem. In other words, the unknowns are not just $\{C_j\}$ and $\phi$, but also $[0, x_d(t))$. Therefore, it is important to specify appropriate boundary conditions that govern the evolution of the domain $[0, x_d(t))$.

First, we specify the initial shape of the domain. The initial conditions are given in the form :

$$\left\{ \begin{array}{l} x_d(0) = x_d^0, \\ C_j(x, 0) = C_j^0(x) \quad \text{in } [0, x_d^0) \quad (j = 1, \dots, M), \end{array} \right. \qquad (2.4)$$

where $x_d^0$ is the initial shape of the pit/crevice as shown in Figure 2.1, and $C_j^0(x)$ ($j = 1, \dots, M$) are functions satisfying

$$C_j^0 \ge 0 \quad (j = 1, \dots, M), \quad z_1 C_1^0 + \dots + z_M C_M^0 = 0. \qquad (2.5)$$

We do not specify the initial value of $\phi$ since $\phi$ is obtained by solving a boundary value problem as we shall see later. Finding appropriate boundary conditions is difficult, which reflects the complex mechanism of corrosion. Later, we state precisely the boundary conditions we choose for the concentrations and the potential at the pit/crevice bottom and the pit/crevice entrance.
In what follows we shall always assume that the solution is completely electrically neutral in the entire pit/crevice solution.

Next we multiply the equations in (2.1a) by $z_j F$ and sum them up for $j = 1, \dots, M$. With equation (2.2), this leads to

$$z_1 F \frac{\partial C_j}{\partial t} + \dots + z_M F \frac{\partial C_M}{\partial t} = F \frac{\partial}{\partial t} (z_1 C_1 + \dots + z_M C_M) = 0.$$

Consequently, we obtain from equation (2.1a),

$$\sum_{j=1}^{M} \frac{\partial}{\partial x} \left( D_j \left( z_j F \frac{\partial C_j}{\partial x} + \frac{z_j^2 F^2}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right) = 0.$$

This equation can be rewritten as

$$\frac{\partial}{\partial x} \left( \sigma \frac{\partial \phi}{\partial x} \right) + \sum_{j=1}^{M} \frac{\partial}{\partial x} \left( z_j F D_j \frac{\partial C_j}{\partial x} \right) = 0, \quad t > 0, 0 < x < x_d(t) \quad \text{where} \quad \sigma := \sum_{j=1}^{M} \frac{z_j^2 F^2}{RT} D_j C_j, \quad (2.3)$$

and $\sigma$ is the electrical conductivity of the ionic solution.

Thus, in order to keep the electroneutrality of the solution, $\phi$ has to satisfy (2.3). In the special case where $C_j$ are spatially uniform (i.e. $\frac{\partial C_j}{\partial x} = 0, j = 1, \dots, M$), equation (2.3) reduces to Laplace's equation $\frac{\partial^2 \phi}{\partial x^2} = 0$, hence $\phi$ is harmonic. However, $\phi$ is not harmonic except in such special cases, and thus in general, the Laplace equation cannot be used (see [13] for a general discussion on this topic).

## 2.2. Initial and boundary conditions

Henceforth, the equation system (2.1)–(2.2) is considered. Needless to say, (2.1)–(2.2) cannot be solved without specifying the boundary conditions as well as the initial condition at $t = 0$. Since the domain $[0, x_d(t))$, which represents the pit/crevice, evolves along with $\{C_j\}$ and $\phi$, this is a free boundary problem. In other words, the unknowns are not just $\{C_j\}$ and $\phi$, but also $[0, x_d(t))$. Therefore, it is important to specify appropriate boundary conditions that govern the evolution of the domain $[0, x_d(t))$.

First, we specify the initial shape of the domain. The initial conditions are given in the form :

$$\left\{ \begin{array}{l} x_d(0) = x_d^0, \\ C_j(x, 0) = C_j^0(x) \quad \text{in } [0, x_d^0) \quad (j = 1, \dots, M), \end{array} \right. \qquad (2.4)$$

where $x_d^0$ is the initial shape of the pit/crevice as shown in Figure 2.1, and $C_j^0(x)$ ($j = 1, \dots, M$) are functions satisfying

$$C_j^0 \ge 0 \quad (j = 1, \dots, M), \quad z_1 C_1^0 + \dots + z_M C_M^0 = 0. \qquad (2.5)$$

We do not specify the initial value of $\phi$ since $\phi$ is obtained by solving a boundary value problem as we shall see later. Finding appropriate boundary conditions is difficult, which reflects the complex mechanism of corrosion. Later, we state precisely the boundary conditions we choose for the concentrations and the potential at the pit/crevice bottom and the pit/crevice entrance.
**Boundary conditions satisfied by {$C_j$}$_{j \in (1,...,M)}$ at the pit bottom $x = x_d(t)$**.

During the development of corrosion pit/crevice, an electrochemical reaction occurs at the pit/crevice bottom. There are two types of electrochemical reaction : oxidation reaction (anodic reaction) which corresponds to the metal dissolution and reduction reaction (cathodic reaction) mainly of the chemical species in the environment such as dissolved O₂ or H⁺. For simplicity, in the present model, we only consider a single electrochemical reaction. More precisely, we neglect the cathodic reaction at the pit/crevice bottom and only consider the anodic reaction of metal dissolution. In our case, the metal to be studied is pure iron. So, this anodic reaction is given by:

$$Fe \underset{k_c}{\overset{k_a}{\rightleftharpoons}} Fe^{2+} + 2e^{-}.$$

To derive boundary conditions at the pit/crevice bottom, we need to distinguish Fe²⁺ and other ions such as Na⁺ and Cl⁻. In what follows, we let $j = 1$ represent Fe²⁺, that is,

$$C_1 := C_{Fe^{2+}}, \quad z_1 := z_{Fe^{2+}} = +2.$$

Then the boundary conditions at $x = x_d(t)$ are given by

$$J_1 = -f_{Fe^{2+}} + \frac{dx_d(t)}{dt} C_1 \quad \text{on } x = x_d(t), \qquad (2.6a)$$

$$J_j = \frac{dx_d(t)}{dt} C_j \quad \text{on } x = x_d(t) \quad (j = 2, \dots, M), \qquad (2.6b)$$

where

$J_j = -D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right)$, $j = 1, \dots, M$, are the ion fluxes, $\frac{dx_d(t)}{dt}$ denotes the velocity of the free boundary $x_d(t)$ at the pit bottom, which we call the *corrosion rate*. It indicates how fast the pit/crevice grows.

The boundary condition (2.6b) is the no-flux boundary condition for the case of a moving boundary. This condition implies that there is no inflow or outflow of ions through the boundary except for iron ions Fe²⁺. On the other hand, the boundary condition (2.6a) implies that there is an inflow of iron ions from the boundary. This inflow $f_{Fe^{2+}}$ denotes the *ferrous ion flux*, that is, the amount of the ferrous ions dissolving from the surface of the iron in the pit/crevice caused by the anodic reaction as mentioned above (per unit time and unit area) minus the amount of the iron ions that is reduced in the pit/crevice solution. For $f_{Fe^{2+}}$, we adopt the following rate expression giving the usual *Butler-Volmer formula* (see [13, p.203-212]) for the oxidation of Fe into Fe²⁺:

$$f_{Fe^{2+}}(\phi, C_{Fe^{2+}}) = k_a \exp\left(\frac{F}{RT}(\phi_m - \phi)\right) - k_c \frac{C_{Fe^{2+}}}{C_{ref}} \exp\left(-\frac{F}{RT}(\phi_m - \phi)\right), \quad (2.7)$$

where $k_a$ and $k_c$ are positive constants that represent the speed of oxidation of iron and that of reduction, respectively (mol/m²s). $C_{ref}$ is a reference concentration which is usually set at $C_{ref} = 1$ mol/L in ideal system, while $\phi_m$ denotes the electric potential applied on the iron steel surface (which is assumed to be constant) and $\phi$, the electric potential of the aqueous solution in contact with the steel surface. In this case, this pit/crevice corrosion phenomenon is considered to be under *potentiostatic conditions*.
Boundary conditions satisfied by {$C_j$}$_{j \in (1,...,M)}$ at the pit bottom $x = x_d(t)$.

During the development of corrosion pit/crevice, an electrochemical reaction occurs at the pit/crevice bottom. There are two types of electrochemical reaction : oxidation reaction (anodic reaction) which corresponds to the metal dissolution and reduction reaction (cathodic reaction) mainly of the chemical species in the environment such as dissolved O₂ or H⁺. For simplicity, in the present model, we only consider a single electrochemical reaction. More precisely, we neglect the cathodic reaction at the pit/crevice bottom and only consider the anodic reaction of metal dissolution. In our case, the metal to be studied is pure iron. So, this anodic reaction is given by:

$$
\mathrm{Fe} \underset{k_c}{\overset{k_a}{\rightleftharpoons}} \mathrm{Fe}^{2+} + 2 \mathrm{e}^{-}.
$$

To derive boundary conditions at the pit/crevice bottom, we need to distinguish Fe²⁺ and other
ions such as Na⁺ and Cl⁻. In what follows, we let j = 1 represent Fe²⁺, that is,

$$
C_1 := C_{\text{Fe}^{2+}}, \quad z_1 := z_{\text{Fe}^{2+}} = +2.
$$

Then the boundary conditions at $x = x_d(t)$ are given by

$$
J_1 = -f_{\text{Fe}^{2+}} + \frac{dx_d(t)}{dt} C_1 \quad \text{on } x = x_d(t), \qquad (2.6a)
$$

$$
\begin{equation}
\begin{cases}
J_j = \frac{dx_d(t)}{dt} C_j & \text{on } x = x_d(t) \quad (j = 2, \dots, M), \\
0 & \text{on } x = x_d(t) \quad (j = 1, \dots, M),
\end{cases}
\tag{2.6b}
\end{equation}
$$

where

$J_j = -D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right)$, $j = 1, \dots, M$, are the ion fluxes,

$\frac{dx_d(t)}{dt}$ denotes the velocity of the free boundary $x_d(t)$ at the pit bottom,

which we call the *corrosion rate*. It indicates how fast the pit/crevice grows.

The boundary condition (2.6b) is the no-flux boundary condition for the case of a moving
boundary. This condition implies that there is no inflow or outflow of ions through the boundary
except for iron ions Fe²⁺. On the other hand, the boundary condition (2.6a) implies that there is
an inflow of iron ions from the boundary. This inflow f<sub>Fe2+</sub> denotes the ferrous ion flux, that is,
the amount of the ferrous ions dissolving from the surface of the iron in the pit/crevice caused by
the anodic reaction as mentioned above (per unit time and unit area) minus the amount of the iron
ions that is reduced in the pit/crevice solution. For f<sub>Fe2+</sub>, we adopt the following rate expression
giving the usual *Butler-Volmer formula* (see [13, p.203-212]) for the oxidation of Fe into Fe²⁺:

$$
f_{\text{Fe}^{2+}}(\phi, C_{\text{Fe}^{2+}}) = k_a \exp\left(\frac{F}{RT}(\phi_m - \phi)\right) - k_c \frac{C_{\text{Fe}^{2+}}}{C_{\text{ref}}} \exp\left(-\frac{F}{RT}(\phi_m - \phi)\right), \quad (2.7)
$$

where $k_a$ and $k_c$ are positive constants that represent the speed of oxidation of iron and that
of reduction, respectively (mol/m²s). $C_{\text{ref}}$ is a reference concentration which is usually set at
$C_{\text{ref}} = 1$ mol/L in ideal system, while $\phi_m$ denotes the electric potential applied on the iron steel
surface (which is assumed to be constant) and $\phi$, the electric potential of the aqueous solution in
contact with the steel surface. In this case, this pit/crevice corrosion phenomenon is considered to
be under *potentiostatic conditions*.
**Boundary conditions satisfied by $\phi$ at the pit/crevice bottom $x = x_d(t)$**.

We can derive the boundary condition for $\phi$ on $x = x_d(t)$ from those for $\{C_j\}$ and the local electroneutrality condition (2.2). Let us multiply (2.6a), (2.6b) by $z_j F$ ($j = 1, \dots, M$) and sum them up, and use (2.2). Then we obtain

$$ \sigma \frac{\partial \phi}{\partial x} + \sum_{j=1}^{M} z_j F D_j \frac{\partial C_j}{\partial x} = z_1 F f_{\text{Fe}^{2+}}(\phi, C_1) \quad \text{on } x = x_d(t), \qquad (2.8) $$

where $\sigma = \sum_{j=1}^{M} \frac{z_j^2 F^2}{RT} D_j C_j$. Here $C_j$, $j = 1, \dots, M$, are regarded as given functions. This is a nonlinear Robin boundary condition for $\phi$.

**Boundary conditions satisfied by $\{C_j\}_{j \in (1, \dots, M)}$ at the entrance $x = 0$ of the pit/crevice.**

For the boundary conditions on $x = 0$ (the entrance of the pit/crevice), we suppose that the solution outside the pit/crevice is well stirred so that the concentration of each ion is always constant on $x = 0$. This implies the Dirichlet boundary conditions on $x = 0$, which are given in the form

$$ C_j = \beta_j \quad \text{at } x = 0, t > 0 \quad (j = 1, 2, \dots, M), \qquad (2.9) $$

where $\beta_j$ ($j = 1, \dots, M$) are constants which satisfy the local electroneutrality condition $\sum_{j=1}^{M} z_j \beta_j = 0$.

**Boundary conditions satisfied by $\phi$ at the entrance $x = 0$ of the pit/crevice.**

Unlike the pit/crevice bottom, where we do not specify the boundary condition for $\phi$, we need to specify the boundary condition for $\phi$ at $x = 0$ in order to determine $\phi$ in $[0, x_d(t))$. The value of the potential at the pit entrance can be complex if there are porous precipitate caps or other deposits, but in the absence of such caps, it gives a good approximation to set the value as follows (see [1]):

$$ \phi = 0 \quad \text{at } x = 0, t > 0. \qquad (2.10) $$

**Boundary condition at the moving interface $x = x_d(t)$**.

The corrosion rate $\frac{dx_d(t)}{dt}$, namely the normal velocity of the boundary $x_d(t)$ in the outward direction, is proportional to the flux $f_{\text{Fe}^{2+}}$ through the following relationship :

$$ \frac{dx_d(t)}{dt} = x'_d(t) = \Omega_{\text{Fe}} f_{\text{Fe}^{2+}} \quad (\text{m/s}), \qquad (2.11) $$

where $\Omega_{\text{Fe}}$ is the molar volume of the solid iron, $\Omega_{\text{Fe}} \approx 7.1 \text{ cm}^3/\text{mol} = 7.1 \cdot 10^{-6} \text{ m}^3/\text{mol}$.

2.3. The full one-dimensional model with Dirichlet boundary conditions at the entrance of the pit/crevice.

In view of Subsections 2.1 and 2.2, the dimensional systems satisfied by the potential $\phi$ and the concentrations $\{C_j\}_{j \in (1, \dots, M)}$ are given by:
**Boundary conditions satisfied by $\phi$ at the pit/crevice bottom $x = x_d(t)$**.

We can derive the boundary condition for $\phi$ on $x = x_d(t)$ from those for $\{C_j\}$ and the local electroneutrality condition (2.2). Let us multiply (2.6a), (2.6b) by $z_j F$ ($j = 1, \dots, M$) and sum them up, and use (2.2). Then we obtain

$$ \sigma \frac{\partial \phi}{\partial x} + \sum_{j=1}^{M} z_j F D_j \frac{\partial C_j}{\partial x} = z_1 F f_{\text{Fe}^{2+}}(\phi, C_1) \quad \text{on } x = x_d(t), \qquad (2.8) $$

where $\sigma = \sum_{j=1}^{M} \frac{z_j^2 F^2}{RT} D_j C_j$. Here $C_j$, $j = 1, \dots, M$, are regarded as given functions. This is a nonlinear Robin boundary condition for $\phi$.

**Boundary conditions satisfied by $\{C_j\}_{j \in (1, \dots, M)}$ at the entrance $x = 0$ of the pit/crevice.**

For the boundary conditions on $x = 0$ (the entrance of the pit/crevice), we suppose that the solution outside the pit/crevice is well stirred so that the concentration of each ion is always constant on $x = 0$. This implies the Dirichlet boundary conditions on $x = 0$, which are given in the form

$$ C_j = \beta_j \quad \text{at } x = 0, t > 0 \quad (j = 1, 2, \dots, M), \qquad (2.9) $$

where $\beta_j$ ($j = 1, \dots, M$) are constants which satisfy the local electroneutrality condition $\sum_{j=1}^{M} z_j \beta_j = 0$.

**Boundary conditions satisfied by $\phi$ at the entrance $x = 0$ of the pit/crevice.**

Unlike the pit/crevice bottom, where we do not specify the boundary condition for $\phi$, we need to specify the boundary condition for $\phi$ at $x = 0$ in order to determine $\phi$ in $[0, x_d(t))$. The value of the potential at the pit entrance can be complex if there are porous precipitate caps or other deposits, but in the absence of such caps, it gives a good approximation to set the value as follows (see [1]):

$$ \phi = 0 \quad \text{at } x = 0, t > 0. \qquad (2.10) $$

**Boundary condition at the moving interface $x = x_d(t)$**.

The corrosion rate $\frac{dx_d(t)}{dt}$, namely the normal velocity of the boundary $x_d(t)$ in the outward direction, is proportional to the flux $f_{\text{Fe}^{2+}}$ through the following relationship :

$$ \frac{dx_d(t)}{dt} = x'_d(t) = \Omega_{\text{Fe}} f_{\text{Fe}^{2+}} \quad (\text{m/s}), \qquad (2.11) $$

where $\Omega_{\text{Fe}}$ is the molar volume of the solid iron, $\Omega_{\text{Fe}} \approx 7.1 \text{ cm}^3/\text{mol} = 7.1 \cdot 10^{-6} \text{ m}^3/\text{mol}$.

2.3. The full one-dimensional model with Dirichlet boundary conditions at the entrance of the pit/crevice.

In view of Subsections 2.1 and 2.2, the dimensional systems satisfied by the potential $\phi$ and the concentrations $\{C_j\}_{j \in (1, \dots, M)}$ are given by:
System satisfied by the potential $\phi$

$$ \left\{ \frac{\partial}{\partial x} \left( \sigma \frac{\partial \phi}{\partial x} \right) + \sum_{j=1}^{M} \frac{\partial}{\partial x} \left( z_j F D_j \frac{\partial C_j}{\partial x} \right) = 0, \quad t > 0, \quad 0 < x < x_d(t), \quad (2.12a) \right. $$

$$ \sigma \frac{\partial \phi}{\partial x} + \sum_{j=1}^{M} z_j F D_j \frac{\partial C_j}{\partial x} = z_1 F f_{\text{Fe}^{2+}}(\phi, C_1) \quad \text{at } x = x_d(t), \quad (2.12b) $$

$$ \phi(0, t) = 0, \quad t > 0, \quad (2.12c) $$

where $\sigma$ is given in (2.3) and $f_{\text{Fe}^{2+}}$ is given in (2.7).

System satisfied by the concentrations $\{C_j\}_{j=1, \dots, M}$

$$ \left\{ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right), \quad t > 0, \quad 0 < x < x_d(t), \quad j = 1, \dots, M, \quad (2.13a) \right. $$

$$ D_1 \left( \frac{\partial C_1}{\partial x} + \frac{z_1 F}{RT} C_1 \frac{\partial \phi}{\partial x} \right) = -x'_d(t) C_1 + f_{\text{Fe}^{2+}}(\phi, C_1) \quad \text{at } x = x_d(t), \quad (2.13b) $$

$$ D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) = -x'_d(t) C_j \quad \text{at } x = x_d(t) \quad (j = 2, \dots, M), \quad (2.13c) $$

$$ C_j(0, t) = \beta_j, \quad t > 0 \quad (j = 1, 2, \dots, M), \quad (2.13d) $$

$$ C_j(x, 0) = C_j^0(x), \quad x \in [0, x_d^0) \quad (2.13e) $$

along with the free boundary condition

$$ x'_{d}(t) = \Omega_{\text{Fe}} f_{\text{Fe}^{2+}} \left( \phi(x_{d}(t), t), C_{1}(x_{d}(t), t) \right) \quad \text{at } x = x_{d}(t), \quad x_{d}(0) = x_{d}^{0}. \quad (2.14) $$

To better understand the behavior of the above mathematical system (2.12)–(2.14), it is convenient to make these equations dimensionless. This latter non-dimensional system is expressed with different symbols: the $\hat{C}_j$’s denote dimensionless concentrations, obtained by normalizing the molar concentrations with a characteristic reference value $C_0$. The independent variables time and space are also normalized by appropriate characteristic values $\tau$ and $L$ respectively, where $\tau = L^2/D_0$ and $D_0$ denotes a characteristic diffusivity, $L$ being a characteristic length. The dimensionless potential is normalized by the ratio $RT/F$. Let us introduce the following changes of variables

$$ x = L\hat{x}, \quad t = \frac{L^2}{D_0}\hat{t}, \quad C_j = C_0\hat{C}_j, \quad D_j = D_0\hat{D}_j, \quad \phi = \frac{RT}{F}\hat{\phi}. \quad (2.15) $$

The system (2.12)–(2.14) then becomes:

System satisfied by the non-dimensional potential $\hat{\phi}$

$$ \left\{ \frac{\partial}{\partial \hat{x}} \left( \hat{\sigma} \frac{\partial \hat{\phi}}{\partial \hat{x}} \right) + \sum_{j=1}^{M} \frac{\partial}{\partial \hat{x}} \left( z_j \hat{D}_j \frac{\partial \hat{C}_j}{\partial \hat{x}} \right) = 0, \quad \hat{t} > 0, \quad 0 < \hat{x} < \hat{x}_d(\hat{t}), \quad (2.16a) \right. $$

$$ \hat{\sigma} \frac{\partial \hat{\phi}}{\partial \hat{x}} + \sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial \hat{C}_j}{\partial \hat{x}} = z_1 \hat{f}_{\text{Fe}^{2+}}(\hat{\phi}, \hat{C}_1) \quad \text{at } \hat{x} = \hat{x}_d(\hat{t}), \quad \hat{C}_1 = \hat{C}_{\text{Fe}^{2+}}, \quad (2.16b) $$

$$ \hat{\phi}(0, \hat{t}) = 0, \quad \hat{t} > 0, \quad (2.16c) $$
System satisfied by the potential $\phi$

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial x} \left( \sigma \frac{\partial \phi}{\partial x} \right) + \sum_{j=1}^{M} \frac{\partial}{\partial x} \left( z_j F D_j \frac{\partial C_j}{\partial x} \right) = 0, \quad t > 0, \ 0 < x < x_d(t), \qquad (2.12a) \\ \displaystyle \sigma \frac{\partial \phi}{\partial x} + \sum_{j=1}^{M} z_j F D_j \frac{\partial C_j}{\partial x} = z_1 F f_{\text{Fe}^{2+}}(\phi, C_1) \quad \text{at } x = x_d(t), \qquad (2.12b) \\ \phi(0, t) = 0, \quad t > 0, \qquad (2.12c) \end{array} \right. $$

where $\sigma$ is given in (2.3) and $f_{\text{Fe}^{2+}}$ is given in (2.7).

System satisfied by the concentrations $\{C_j\}_{j=1, \dots, M}$

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right), \quad t > 0, \ 0 < x < x_d(t), \quad j = 1, \dots, M, \qquad (2.13a) \\ D_1 \left( \frac{\partial C_1}{\partial x} + \frac{z_1 F}{RT} C_1 \frac{\partial \phi}{\partial x} \right) = -x'_d(t) C_1 + f_{\text{Fe}^{2+}}(\phi, C_1) \quad \text{at } x = x_d(t), \qquad (2.13b) \\ D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) = -x'_d(t) C_j \quad \text{at } x = x_d(t) \quad (j = 2, \dots, M), \qquad (2.13c) \\ C_j(0, t) = \beta_j, \quad t > 0 \quad (j = 1, 2, \dots, M), \qquad (2.13d) \\ C_j(x, 0) = C_j^0(x), \quad x \in [0, x_d^0] \qquad (2.13e) \end{array} \right. $$

along with the free boundary condition

$$ x'_{d}(t) = \Omega_{\text{Fe}} f_{\text{Fe}^{2+}} \left( \phi(x_d(t), t), C_1(x_d(t), t) \right) \quad \text{at } x = x_d(t), \quad x_d(0) = x_d^0. \qquad (2.14) $$

To better understand the behavior of the above mathematical system (2.12)–(2.14), it is convenient to make these equations dimensionless. This latter non-dimensional system is expressed with different symbols: the $\hat{C}_j$’s denote dimensionless concentrations, obtained by normalizing the molar concentrations with a characteristic reference value $C_0$. The independent variables time and space are also normalized by appropriate characteristic values $\tau$ and $L$ respectively, where $\tau = L^2/D_0$ and $D_0$ denotes a characteristic diffusivity, $L$ being a characteristic length. The dimensionless potential is normalized by the ratio $RT/F$. Let us introduce the following changes of variables

$$ x = L\hat{x}, \quad t = \frac{L^2}{D_0}\hat{t}, \quad C_j = C_0\hat{C}_j, \quad D_j = D_0\hat{D}_j, \quad \phi = \frac{RT}{F}\hat{\phi}. \qquad (2.15) $$

The system (2.12)–(2.14) then becomes:

System satisfied by the non-dimensional potential $\hat{\phi}$

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial \hat{x}} \left( \hat{\sigma} \frac{\partial \hat{\phi}}{\partial \hat{x}} \right) + \sum_{j=1}^{M} \frac{\partial}{\partial \hat{x}} \left( z_j \hat{D}_j \frac{\partial \hat{C}_j}{\partial \hat{x}} \right) = 0, \quad \hat{t} > 0, \ 0 < \hat{x} < \hat{x}_d(\hat{t}), \qquad (2.16a) \\ \displaystyle \hat{\sigma} \frac{\partial \hat{\phi}}{\partial \hat{x}} + \sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial \hat{C}_j}{\partial \hat{x}} = z_1 \hat{f}_{\text{Fe}^{2+}}(\hat{\phi}, \hat{C}_1) \quad \text{at } \hat{x} = \hat{x}_d(\hat{t}), \quad \hat{C}_1 = \hat{C}_{\text{Fe}^{2+}}, \qquad (2.16b) \\ \displaystyle \hat{\phi}(0, \hat{t}) = 0, \quad \hat{t} > 0, \qquad (2.16c) \end{array} \right. $$
where $\hat{\sigma} := \sum_{j=1}^{M} z_j^2 \hat{D}_j \hat{C}_j$ and the non-dimensional form of $f_{\text{Fe}^{2+}}$ denoted by $\hat{f}_{\text{Fe}^{2+}}$ is given by

$$ \hat{f}_{\text{Fe}^{2+}}(\hat{\phi}, \hat{C}_1) = \hat{k}_a \exp(\hat{\phi}_m - \hat{\phi}) - \hat{k}_c \hat{C}_1 \exp(-(\hat{\phi}_m - \hat{\phi})), \quad (2.17) $$

where

$$ \hat{\phi}_m = \frac{F}{RT}\phi_m, \quad \hat{k}_a = \frac{L}{D_0 C_0} k_a \quad \text{and} \quad \hat{k}_c = \frac{L}{D_0 C_{\text{ref}}} k_c. \quad (2.18) $$

System satisfied by the non-dimensional concentrations $\{\hat{C}_j\}_{j=1,\dots,M}$

$$ \left\{ \frac{\partial \hat{C}_j}{\partial \hat{t}} = \frac{\partial}{\partial \hat{x}} \left( \hat{D}_j \left( \frac{\partial \hat{C}_j}{\partial \hat{x}} + z_j \hat{C}_j \frac{\partial \hat{\phi}}{\partial \hat{x}} \right) \right), \quad \hat{t} > 0, 0 < \hat{x} < \hat{x}_d(\hat{t}), \quad j = 1, \dots, M, \quad (2.19a) $$

$$ \hat{D}_1 \left( \frac{\partial \hat{C}_1}{\partial \hat{x}} + z_1 \hat{C}_1 \frac{\partial \hat{\phi}}{\partial \hat{x}} \right) = -\hat{x}'_d(\hat{t}) \hat{C}_1 + \hat{f}_{\text{Fe}^{2+}}(\hat{\phi}, \hat{C}_1) \quad \text{at} \quad \hat{x} = \hat{x}_d(\hat{t}), \quad (2.19b) $$

$$ \hat{D}_j \left( \frac{\partial \hat{C}_j}{\partial \hat{x}} + z_j \hat{C}_j \frac{\partial \hat{\phi}}{\partial \hat{x}} \right) = -\hat{x}'_d(\hat{t}) \hat{C}_j \quad \text{at} \quad \hat{x} = \hat{x}_d(\hat{t}) \quad (j = 2, \dots, M), \quad (2.19c) $$

$$ \hat{C}_j(0, \hat{t}) = \hat{\beta}_j, \quad \hat{t} > 0 \quad (j = 1, 2, \dots, M) \quad \text{where} \quad \hat{\beta}_j = \frac{\beta_j}{C_0}, \quad (2.19d) $$

$$ C_j(x, 0) = C_j^0(x), \quad x \in [0, x_d^0) \quad (2.19e) $$

along with the non-dimensional free boundary condition

$$ \hat{x}'_d(\hat{t}) = \hat{\Omega}_{\text{Fe}} \hat{f}_{\text{Fe}^{2+}} \left( \hat{\phi}(\hat{x}_d(\hat{t}), \hat{t}), \hat{C}_1(\hat{x}_d(\hat{t}), \hat{t}) \right) \quad \text{at} \quad \hat{x} = \hat{x}_d(\hat{t}) \quad \text{where} \quad \hat{\Omega}_{\text{Fe}} = C_0 \Omega_{\text{Fe}}, \quad x_d(0) = x_d^0. \quad (2.20) $$

Later, the analysis will be done for the non-dimensionalized system (2.16)–(2.19),(2.20), but for notational simplicity we will drop the symbol $\hat{\cdot}$ from variables $\hat{C}_j, \hat{\phi}, \hat{x}, \hat{t}$. Henceforth, the following non-dimensional system is considered from now on:

$$ \left\{ \frac{\partial}{\partial x} \left( \hat{\sigma} \frac{\partial \hat{\phi}}{\partial x} \right) + \sum_{j=1}^{M} \frac{\partial}{\partial x} \left( z_j \hat{D}_j \frac{\partial C_j}{\partial x} \right) = 0, \quad t > 0, 0 < x < x_d(t), \quad (2.21a) $$

$$ \hat{\sigma} \frac{\partial \hat{\phi}}{\partial x} + \sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial C_j}{\partial x} = z_1 \hat{f}_{\text{Fe}^{2+}}(\hat{\phi}, C_1) \quad \text{at} \quad x = x_d(t), \quad (2.21b) $$

$$ \phi(0, t) = 0, \quad t > 0, \quad (2.21c) $$

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( \hat{D}_j \left( \frac{\partial C_j}{\partial x} + z_j C_j \frac{\partial \hat{\phi}}{\partial x} \right) \right), \quad t > 0, 0 < x < x_d(t), \quad (j = 1, \dots, M) \quad (2.22a) $$

$$ \hat{D}_1 \frac{\partial C_1}{\partial x} + z_1 \hat{D}_1 C_1 \frac{\partial \hat{\phi}}{\partial x} = -x'_d(t) C_1 + \hat{f}_{\text{Fe}^{2+}}(\hat{\phi}, C_1) \quad \text{at} \quad x = x_d(t), \quad (2.22b) $$

$$ \hat{D}_j \frac{\partial C_j}{\partial x} + z_j \hat{D}_j C_j \frac{\partial \hat{\phi}}{\partial x} = -x'_d(t) C_j \quad \text{at} \quad x = x_d(t) \quad (j = 2, \dots, M), \quad (2.22c) $$

$$ C_j(0, t) = \hat{\beta}_j, \quad t > 0 \quad (j = 1, \dots, M), \quad (2.22d) $$

$$ C_j(x, 0) = C_j^0(x), \quad x \in [0, x_d^0), \quad (2.22e) $$
where $\hat{\sigma} := \sum_{j=1}^{M} z_j^2 \hat{D}_j \hat{C}_j$ and the non-dimensional form of $f_{\text{Fe}^{2+}}$ denoted by $\hat{f}_{\text{Fe}^{2+}}$ is given by

$$ \hat{f}_{\text{Fe}^{2+}}(\hat{\phi}, \hat{C}_1) = \hat{k}_a \exp(\hat{\phi}_m - \hat{\phi}) - \hat{k}_c \hat{C}_1 \exp(-(\hat{\phi}_m - \hat{\phi})), \quad (2.17) $$

where

$$ \hat{\phi}_m = \frac{F}{RT}\phi_m, \quad \hat{k}_a = \frac{L}{D_0 C_0} k_a \quad \text{and} \quad \hat{k}_c = \frac{L}{D_0 C_{\text{ref}}} k_c. \quad (2.18) $$

System satisfied by the non-dimensional concentrations $\{\hat{C}_j\}_{j=1,\dots,M}$

$$ \left\{ \frac{\partial \hat{C}_j}{\partial \hat{t}} = \frac{\partial}{\partial \hat{x}} \left( \hat{D}_j \left( \frac{\partial \hat{C}_j}{\partial \hat{x}} + z_j \hat{C}_j \frac{\partial \hat{\phi}}{\partial \hat{x}} \right) \right), \quad \hat{t} > 0, 0 < \hat{x} < \hat{x}_d(\hat{t}), \quad j = 1, \dots, M, \quad (2.19a) $$

$$ \hat{D}_1 \left( \frac{\partial \hat{C}_1}{\partial \hat{x}} + z_1 \hat{C}_1 \frac{\partial \hat{\phi}}{\partial \hat{x}} \right) = -\hat{x}'_d(\hat{t}) \hat{C}_1 + \hat{f}_{\text{Fe}^{2+}}(\hat{\phi}, \hat{C}_1) \quad \text{at} \quad \hat{x} = \hat{x}_d(\hat{t}), \quad (2.19b) $$

$$ \hat{D}_j \left( \frac{\partial \hat{C}_j}{\partial \hat{x}} + z_j \hat{C}_j \frac{\partial \hat{\phi}}{\partial \hat{x}} \right) = -\hat{x}'_d(\hat{t}) \hat{C}_j \quad \text{at} \quad \hat{x} = \hat{x}_d(\hat{t}) \quad (j = 2, \dots, M), \quad (2.19c) $$

$$ \hat{C}_j(0, \hat{t}) = \hat{\beta}_j, \quad \hat{t} > 0 \quad (j = 1, 2, \dots, M) \quad \text{where} \quad \hat{\beta}_j = \frac{\beta_j}{C_0}, \quad (2.19d) $$

$$ C_j(x, 0) = C_j^0(x), \quad x \in [0, x_d^0) \quad (2.19e) $$

along with the non-dimensional free boundary condition

$$ \hat{x}'_d(\hat{t}) = \hat{\Omega}_{\text{Fe}} \hat{f}_{\text{Fe}^{2+}} \left( \hat{\phi}(\hat{x}_d(\hat{t}), \hat{t}), \hat{C}_1(\hat{x}_d(\hat{t}), \hat{t}) \right) \quad \text{at} \quad \hat{x} = \hat{x}_d(\hat{t}) \quad \text{where} \quad \hat{\Omega}_{\text{Fe}} = C_0 \Omega_{\text{Fe}}, \quad x_d(0) = x_d^0. \quad (2.20) $$

Later, the analysis will be done for the non-dimensionalized system (2.16)–(2.19), (2.20), but for notational simplicity we will drop the symbol $\hat{\cdot}$ from variables $\hat{C}_j, \hat{\phi}, \hat{x}, \hat{t}$. Henceforth, the following non-dimensional system is considered from now on:

$$ \left\{ \frac{\partial}{\partial x} \left( \hat{\sigma} \frac{\partial \hat{\phi}}{\partial x} \right) + \sum_{j=1}^{M} \frac{\partial}{\partial x} \left( z_j \hat{D}_j \frac{\partial C_j}{\partial x} \right) = 0, \quad t > 0, 0 < x < x_d(t), \quad (2.21a) $$

$$ \hat{\sigma} \frac{\partial \hat{\phi}}{\partial x} + \sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial C_j}{\partial x} = z_1 \hat{f}_{\text{Fe}^{2+}}(\hat{\phi}, C_1) \quad \text{at} \quad x = x_d(t), \quad (2.21b) $$

$$ \phi(0, t) = 0, \quad t > 0, \quad (2.21c) $$

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( \hat{D}_j \left( \frac{\partial C_j}{\partial x} + z_j C_j \frac{\partial \hat{\phi}}{\partial x} \right) \right), \quad t > 0, \quad 0 < x < x_d(t), \quad (j = 1, \dots, M) \quad (2.22a) $$

$$ \hat{D}_1 \frac{\partial C_1}{\partial x} + z_1 \hat{D}_1 C_1 \frac{\partial \hat{\phi}}{\partial x} = -x'_d(t) C_1 + \hat{f}_{\text{Fe}^{2+}}(\hat{\phi}, C_1) \quad \text{at} \quad x = x_d(t), \quad (2.22b) $$

$$ \hat{D}_j \frac{\partial C_j}{\partial x} + z_j \hat{D}_j C_j \frac{\partial \hat{\phi}}{\partial x} = -x'_d(t) C_j \quad \text{at} \quad x = x_d(t) \quad (j = 2, \dots, M), \quad (2.22c) $$

$$ C_j(0, t) = \hat{\beta}_j, \quad t > 0 \quad (j = 1, \dots, M), \quad (2.22d) $$

$$ C_j(x, 0) = C_j^0(x), \quad x \in [0, x_d^0), \quad (2.22e) $$
along with the free boundary condition

$$x'_{d}(t) = \hat{\Omega}_{\text{Fe}} \hat{f}_{\text{Fe}^{2+}} \left( \phi(x_{d}(t), t), C_{1}(x_{d}(t), t) \right) \quad \text{at } x = x_{d}(t), \quad x_{d}(0) = x_{d}^{0}, \qquad (2.23)$$

where

$$\hat{f}_{\text{Fe}^{2+}}(\phi, C_1) = \hat{k}_a \exp(\hat{\phi}_m - \phi) - \hat{k}_c C_1 \exp(-(\hat{\phi}_m - \phi)). \qquad (2.24)$$

### 3. Numerical approximation of the one-dimensional model

#### 3.1. Introduction

There is a large choice of methods for numerically solving one-dimensional problems which involve a moving interface [14, 15, 16, 17]. One of the popular techniques used is the Arbitrary Lagrangian Eulerian (ALE) method (see [18, 19]). To solve the non-dimensionalized system (2.21)–(2.23), a special choice of an ALE method is applied. In the literature, this method is usually referred to as the Variable Space Grid (VSG) method [16, 20]. Considering the ALE formulation of the diffusion-convection system (2.21)–(2.23), the used numerical method is based on a decoupling of the displacement of the free boundary from the numerical solution of the diffusion-convection equation. The motion of the free boundary is numerically solved by an explicit Euler scheme whereas an implicit finite difference scheme is used for the discretization of the ALE formulation of the diffusion-convection equation. Finally, a Newton method is employed to solve this full discretized implicit system (see [12] for complete details). The ALE formulation of the diffusion-convection system (2.21)–(2.23) and the numerical ALE scheme are detailed in Appendix A and Appendix B respectively.

Recall that the system (2.21) for the potential has been obtained by combining the electroneutrality condition (2.2) with the use of the diffusion-convection system (2.22),(2.23) for concentrations (see Section 2.). Thus, the potential system (2.21) can be replaced by the electroneutrality condition:

$$\left\{ \begin{array}{l} \displaystyle \sum_{j=1}^{M} z_j C_j(x,t) = 0, \quad t > 0, \quad 0 \le x \le x_d(t), \qquad (3.1a) \\ \phi(0,t) = 0, \quad t > 0. \qquad (3.1b) \end{array} \right.$$

Hence solving the coupled system (2.21),(2.22),(2.23) is mathematically equivalent to solving the coupled system (3.1),(2.22),(2.23), but from the numerical point of view, the latter is easier to handle. Therefore, in what follows, we shall derive the numerical scheme to solve the coupled system (3.1),(2.22),(2.23) instead of (2.21),(2.22),(2.23). However, the system (2.21) is used at initial time ($t=0$) to determine the initial profile of potential. The next subsection explains this computation.

#### 3.2. Procedure to determine the initial profile of the potential in the pit/crevice

The procedure for solving the system (2.21) to determine the initial profile of the potential $\phi$ is described below. At $t = 0$, the only unknown in (2.21) is the potential $\phi$ since $C_1, C_2, \dots, C_M, x_d$ are given functions at this stage. In fact, for all $t \ge 0$, the potential $\phi$ can be determined as a function of $C_1, C_2, \dots, C_M$ and $x_d$. Indeed, integrating (2.21a) by $x$ shows that

$$\hat{\sigma} \frac{\partial \phi}{\partial x} + \sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial C_j}{\partial x} \qquad (3.2)$$
along with the free boundary condition

$$x'_{d}(t) = \hat{\Omega}_{\text{Fe}} \hat{f}_{\text{Fe}^{2+}} \left( \phi(x_{d}(t), t), C_{1}(x_{d}(t), t) \right) \quad \text{at } x = x_{d}(t), \quad x_{d}(0) = x_{d}^{0}, \qquad (2.23)$$

where

$$\hat{f}_{\text{Fe}^{2+}}(\phi, C_1) = \hat{k}_a \exp(\hat{\phi}_m - \phi) - \hat{k}_c C_1 \exp(-(\hat{\phi}_m - \phi)). \qquad (2.24)$$

## 3. Numerical approximation of the one-dimensional model

### 3.1. Introduction

There is a large choice of methods for numerically solving one-dimensional problems which involve a moving interface [14, 15, 16, 17]. One of the popular techniques used is the Arbitrary Lagrangian Eulerian (ALE) method (see [18, 19]). To solve the non-dimensionalized system (2.21)–(2.23), a special choice of an ALE method is applied. In the literature, this method is usually referred to as the Variable Space Grid (VSG) method [16, 20]. Considering the ALE formulation of the diffusion-convection system (2.21)–(2.23), the used numerical method is based on a decoupling of the displacement of the free boundary from the numerical solution of the diffusion-convection equation. The motion of the free boundary is numerically solved by an explicit Euler scheme whereas an implicit finite difference scheme is used for the discretization of the ALE formulation of the diffusion-convection equation. Finally, a Newton method is employed to solve this full discretized implicit system (see [12] for complete details). The ALE formulation of the diffusion-convection system (2.21)–(2.23) and the numerical ALE scheme are detailed in Appendix A and Appendix B respectively.

Recall that the system (2.21) for the potential has been obtained by combining the electroneutrality condition (2.2) with the use of the diffusion-convection system (2.22),(2.23) for concentrations (see Section 2.). Thus, the potential system (2.21) can be replaced by the electroneutrality condition:

$$\left\{ \begin{array}{l} \displaystyle \sum_{j=1}^{M} z_j C_j(x, t) = 0, \quad t > 0, \quad 0 \le x \le x_d(t), \\ \phi(0, t) = 0, \quad t > 0. \end{array} \right. \qquad (3.1a) \\ (3.1b)$$

Hence solving the coupled system (2.21),(2.22),(2.23) is mathematically equivalent to solving the coupled system (3.1),(2.22),(2.23), but from the numerical point of view, the latter is easier to handle. Therefore, in what follows, we shall derive the numerical scheme to solve the coupled system (3.1),(2.22),(2.23) instead of (2.21),(2.22),(2.23). However, the system (2.21) is used at initial time (t=0) to determine the initial profile of potential. The next subsection explains this computation.

### 3.2. Procedure to determine the initial profile of the potential in the pit/crevice

The procedure for solving the system (2.21) to determine the initial profile of the potential $\phi$ is described below. At $t = 0$, the only unknown in (2.21) is the potential $\phi$ since $C_1, C_2, \dots, C_M, x_d$ are given functions at this stage. In fact, for all $t \ge 0$, the potential $\phi$ can be determined as a function of $C_1, C_2, \dots, C_M$ and $x_d$. Indeed, integrating (2.21a) by $x$ shows that

$$\hat{\sigma} \frac{\partial \phi}{\partial x} + \sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial C_j}{\partial x} \qquad (3.2)$$
is independent of $x$, hence, by (2.21b), the following equation is obtained

$$ \hat{\sigma} \frac{\partial \phi}{\partial x} + \sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial C_j}{\partial x} = z_1 \hat{f}_{\text{Fe}^{2+}}(P(t), C_1(x_d(t), t)) \quad \text{for } 0 \le x \le x_d(t), t > 0 \quad (3.3) $$

where

$$ P(t) := \phi(x_d(t), t). \quad (3.4) $$

Then (3.3) can be rewritten as

$$ \frac{\partial \phi}{\partial x} = - \frac{\sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial C_j}{\partial x}}{\sum_{j=1}^{M} z_j^2 \hat{D}_j C_j} + \frac{z_1 \hat{f}_{\text{Fe}^{2+}}(P(t), C_1(x_d(t), t))}{\sum_{j=1}^{M} z_j^2 \hat{D}_j C_j} \quad (3.5) $$

Integrating this equality from $x = 0$ to $x = x_d(t)$ and using (2.21c), give the following equation

$$ P(t) = -A(t) + 2 \hat{f}_{\text{Fe}^{2+}}(P(t), C_1(x_d(t), t)) B(t), \quad (3.6) $$

where

$$ A(t) = \int_{0}^{x_d(t)} \frac{\sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial C_j}{\partial x}}{\sum_{j=1}^{M} z_j^2 \hat{D}_j C_j} dx, \quad B(t) = \int_{0}^{x_d(t)} \frac{dx}{\sum_{j=1}^{M} z_j^2 \hat{D}_j C_j} \left( = \int_{0}^{x_d(t)} \frac{dx}{\hat{\sigma}(x, t)} \right). \quad (3.7) $$

Here $A(t)$, $B(t)$, $C_1(x_d(t), t)$ are treated as given functions of the known concentrations and the pit/crevice bottom position.

For fixed $t \ge 0$, the potential at the pit/crevice bottom, $P(t) = \phi(x_d(t), t)$, is determined by solving the nonlinear equation (3.6). Once the value of $P(t)$ is obtained (from the given functions $\{C_j\}_{j \in \{1, \dots, M\}}$ and $x_d$), the derivative $\frac{\partial \phi}{\partial x}$ is determined on the entire interval $0 \le x \le x_d(t)$ from (3.5). Thus, the whole initial potential $\phi^0 = \phi(\cdot, 0)$ in $[0, x_d^0]$ at $t = 0$ can be computed from the initial data $x_d^0$ and $C_j^0$, ($j = 1, \dots, M$).

# 4. Numerical simulations of the unidimensional pit/crevice

## 4.1. Comparison of the numerical model with a steady-state model from the literature

To examine the validity of the above mathematical model, we compare the numerical simulations of our model with those of a similar one-dimensional model discussed in the literature. Sharland and Tasker [21, 7, 6] have developed a predictive mathematical model incorporating the electrochemical, chemical and ionic migration processes. Their model predicts the steady-state chemistry solution and the electrode potential within a corroding crevice or pit.

In the case of propagation of a crevice-type corrosion with passive walls, Sharland et al. [7] solved their model using a number of mathematical approximations based on physical observations. In their model, they considered that the cavity propagation was slow compared to the ionic transport, so that the moving boundary effects was ignored. In other words, the geometry of the crevice was
is independent of $x$, hence, by (2.21b), the following equation is obtained

$$ \hat{\sigma} \frac{\partial \phi}{\partial x} + \sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial C_j}{\partial x} = z_1 \hat{f}_{\text{Fe}^{2+}}(P(t), C_1(x_d(t), t)) \quad \text{for } 0 \le x \le x_d(t), t > 0 \quad (3.3) $$

where

$$ P(t) := \phi(x_d(t), t). \quad (3.4) $$

Then (3.3) can be rewritten as

$$ \frac{\partial \phi}{\partial x} = - \frac{\sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial C_j}{\partial x}}{\sum_{j=1}^{M} z_j^2 \hat{D}_j C_j} + \frac{z_1 \hat{f}_{\text{Fe}^{2+}}(P(t), C_1(x_d(t), t))}{\sum_{j=1}^{M} z_j^2 \hat{D}_j C_j} \quad (3.5) $$

Integrating this equality from $x = 0$ to $x = x_d(t)$ and using (2.21c), give the following equation

$$ P(t) = -A(t) + 2 \hat{f}_{\text{Fe}^{2+}}(P(t), C_1(x_d(t), t)) B(t), \quad (3.6) $$

where

$$ A(t) = \int_{0}^{x_d(t)} \frac{\sum_{j=1}^{M} z_j \hat{D}_j \frac{\partial C_j}{\partial x}}{\sum_{j=1}^{M} z_j^2 \hat{D}_j C_j} dx, \quad B(t) = \int_{0}^{x_d(t)} \frac{dx}{\sum_{j=1}^{M} z_j^2 \hat{D}_j C_j} \left( = \int_{0}^{x_d(t)} \frac{dx}{\hat{\sigma}(x, t)} \right). \quad (3.7) $$

Here $A(t)$, $B(t)$, $C_1(x_d(t), t)$ are treated as given functions of the known concentrations and the pit/crevice bottom position.

For fixed $t \ge 0$, the potential at the pit/crevice bottom, $P(t) = \phi(x_d(t), t)$, is determined by solving the nonlinear equation (3.6). Once the value of $P(t)$ is obtained (from the given functions $\{C_j\}_{j \in \{1, \dots, M\}}$ and $x_d$), the derivative $\frac{\partial \phi}{\partial x}$ is determined on the entire interval $0 \le x \le x_d(t)$ from (3.5). Thus, the whole initial potential $\phi^0 = \phi(\cdot, 0)$ in $[0, x_d^0]$ at $t = 0$ can be computed from the initial data $x_d^0$ and $C_j^0$, ($j = 1, \dots, M$).

# 4. Numerical simulations of the unidimensional pit/crevice

## 4.1. Comparison of the numerical model with a steady-state model from the literature

To examine the validity of the above mathematical model, we compare the numerical simulations of our model with those of a similar one-dimensional model discussed in the literature. Sharland and Tasker [21, 7, 6] have developed a predictive mathematical model incorporating the electrochemical, chemical and ionic migration processes. Their model predicts the steady-state chemistry solution and the electrode potential within a corroding crevice or pit.

In the case of propagation of a crevice-type corrosion with passive walls, Sharland et al. [7] solved their model using a number of mathematical approximations based on physical observations. In their model, they considered that the cavity propagation was slow compared to the ionic transport, so that the moving boundary effects was ignored. In other words, the geometry of the crevice was
kept constant, which was a very rough approximation since a real crevice never reaches a steady state because of its shape evolution [6, 7, 22]. The same simulation was reproduced by Mousson et al. [22] using a commercial code.

To be as close as possible to these studies, the following chemical reaction is also considered in the equation system (2.22):

$$ \mathrm{Fe}^{2+} + \mathrm{H}_2\mathrm{O} \underset{\mathrm{k}_b}{\overset{\mathrm{k}_f}{\rightleftharpoons}} \mathrm{Fe}(\mathrm{OH})^{+} + \mathrm{H}^{+}. \quad (4.1) $$

The water autoionization reaction considered in Sharland's model was neglected here to simplify the chemical system. The water autoionization reaction produces $H^+$ and $OH^-$ species. However, we conjecture that the fact that it is not considered in our model has negligible effects on the concentration of the species in the crevice solution. Indeed, the concentration of $H^+$ is mainly controlled by the hydrolysis reaction of iron cations (4.1) which has a higher thermodynamic constant [6] Moreover, the concentration of $OH^-$ species by the water autoionization reaction would mainly affect the $Cl^-$ concentration profile through the electroneutrality condition. However, the $OH^-$ concentration would be extremely low because of the very acidic environment. In consequence, the predicted $Cl^-$ concentration which reaches very high values (> 10 M/L) in our model would only be very slightly decreased by the consideration of the autoionization reaction. Thus, it appears acceptable to compare the results of our simplified model with Sharland's one.

Therefore, we consider a 1D model with five ions including diffusion, migration and reaction terms. In this case $M = 5$ and the species taken into account are $Fe^{2+}$, $Na^+$, $Cl^-$, $H^+$ and $Fe(OH)^+$. A reaction term coming from the chemical reaction (4.1) must be added in the PDE system (2.21)–(2.23). In Appendix C, we describe the reaction term and the resulting full diffusion-migration-reaction system for 5 ions. We also provide the numerical values of the kinetic parameters which are involved. We refer to [12] for a more detailed description of the model with reaction terms.

Computations are performed with the following choice of the reference parameters:

$$ D_0 = 1 \cdot 10^{-9} \, m^2 \cdot s^{-1}, \quad L = 10^{-3} \, m, \quad C_0 = 10^{-3} \, \text{mol/L}. $$

The diffusion coefficients of concerned species are given by (see [23]):

$$ \begin{gathered} D_{\mathrm{Fe}^{2+}} = 7.2 \cdot 10^{-10} \, m^2 \cdot s^{-1}, \quad D_{\mathrm{Cl}^{-}} = 2 \cdot 10^{-9} \, m^2 \cdot s^{-1}, \quad D_{\mathrm{Na}^{+}} = 1.3 \cdot 10^{-9} \, m^2 \cdot s^{-1}, \\ D_{\mathrm{H}^{+}} = 9.3 \cdot 10^{-9} \, m^2 \cdot s^{-1} \quad \text{and} \quad D_{\mathrm{Fe}(\mathrm{OH})^{+}} = 0.8 \cdot 10^{-9} \, m^2 \cdot s^{-1}. \end{gathered} $$

The initial data used in computation are given by:

$$ \begin{gathered} x_d^0 = 2 \, \text{mm}, \quad C_{\mathrm{Fe}^{2+}}^0(x) = 9.8 \cdot 10^{-7} \, \text{mol/L}, \quad C_{\mathrm{H}^+}^0(x) = C_{\mathrm{Fe}(\mathrm{OH})^+}^0(x) = 2.2 \cdot 10^{-8} \, \text{mol/L}, \\ C_{\mathrm{Na}^+}^0(x) = 10^{-3} \, \text{mol/L} \quad \text{and} \quad C_{\mathrm{Cl}^-}^0(x) \approx 10^{-3} \, \text{mol/L}, \quad \text{for} \ 0 < x < x_d^0, \end{gathered} $$

where $x_d^0$ is the initial depth of the crevice and $C_{\mathrm{Cl}^-}^0$ is given to satisfy the local electroneutrality condition of the solution.

Figure 4.1 shows the solution chemistry obtained by Sharland's steady-state model where the moving boundary is ignored and the pit/crevice bottom position is kept fixed and constant. Figure 4.2 shows the concentration profiles of the ions and the crevice depth calculated by our model at two different times of propagation for the same initial conditions as the ones used in Figure 4.1 from our model.
kept constant, which was a very rough approximation since a real crevice never reaches a steady state because of its shape evolution [6, 7, 22]. The same simulation was reproduced by Mousson et al. [22] using a commercial code.

To be as close as possible to these studies, the following chemical reaction is also considered in the equation system (2.22):

$$ \mathrm{Fe}^{2+} + \mathrm{H}_2\mathrm{O} \underset{\mathrm{k}_b}{\overset{\mathrm{k}_f}{\rightleftharpoons}} \mathrm{Fe}(\mathrm{OH})^{+} + \mathrm{H}^{+}. \quad (4.1) $$

The water autoionization reaction considered in Sharland's model was neglected here to simplify the chemical system. The water autoionization reaction produces $H^+$ and $OH^-$ species. However, we conjecture that the fact that it is not considered in our model has negligible effects on the concentration of the species in the crevice solution. Indeed, the concentration of $H^+$ is mainly controlled by the hydrolysis reaction of iron cations (4.1) which has a higher thermodynamic constant [6] Moreover, the concentration of $OH^-$ species by the water autoionization reaction would mainly affect the $Cl^-$ concentration profile through the electroneutrality condition. However, the $OH^-$ concentration would be extremely low because of the very acidic environment. In consequence, the predicted $Cl^-$ concentration which reaches very high values (> 10 M/L) in our model would only be very slightly decreased by the consideration of the autoionization reaction. Thus, it appears acceptable to compare the results of our simplified model with Sharland's one.

Therefore, we consider a 1D model with five ions including diffusion, migration and reaction terms. In this case $M = 5$ and the species taken into account are $Fe^{2+}$, $Na^+$, $Cl^-$, $H^+$ and $Fe(OH)^+$. A reaction term coming from the chemical reaction (4.1) must be added in the PDE system (2.21)–(2.23). In Appendix C, we describe the reaction term and the resulting full diffusion-migration-reaction system for 5 ions. We also provide the numerical values of the kinetic parameters which are involved. We refer to [12] for a more detailed description of the model with reaction terms.

Computations are performed with the following choice of the reference parameters:

$$ D_0 = 1 \cdot 10^{-9} \, \mathrm{m}^2 \cdot \mathrm{s}^{-1}, \quad L = 10^{-3} \, \mathrm{m}, \quad C_0 = 10^{-3} \, \mathrm{mol/L}. $$

The diffusion coefficients of concerned species are given by (see [23]):

$$ \begin{gathered} D_{\mathrm{Fe}^{2+}} = 7.2 \cdot 10^{-10} \, \mathrm{m}^2 \cdot \mathrm{s}^{-1}, \quad D_{\mathrm{Cl}^{-}} = 2 \cdot 10^{-9} \, \mathrm{m}^2 \cdot \mathrm{s}^{-1}, \quad D_{\mathrm{Na}^{+}} = 1.3 \cdot 10^{-9} \, \mathrm{m}^2 \cdot \mathrm{s}^{-1}, \\ D_{\mathrm{H}^{+}} = 9.3 \cdot 10^{-9} \, \mathrm{m}^2 \cdot \mathrm{s}^{-1} \quad \text{and} \quad D_{\mathrm{Fe}(\mathrm{OH})^{+}} = 0.8 \cdot 10^{-9} \, \mathrm{m}^2 \cdot \mathrm{s}^{-1}. \end{gathered} $$

The initial data used in computation are given by:

$$ \begin{gathered} x_d^0 = 2 \, \mathrm{mm}, \quad C_{\mathrm{Fe}^{2+}}^0(x) = 9.8 \cdot 10^{-7} \, \mathrm{mol/L}, \quad C_{\mathrm{H}^+}^0(x) = C_{\mathrm{Fe}(\mathrm{OH})^+}^0(x) = 2.2 \cdot 10^{-8} \, \mathrm{mol/L}, \\ C_{\mathrm{Na}^+}^0(x) = 10^{-3} \, \mathrm{mol/L} \quad \text{and} \quad C_{\mathrm{Cl}^-}^0(x) \approx 10^{-3} \, \mathrm{mol/L}, \quad \text{for} \ 0 < x < x_d^0, \end{gathered} $$

where $x_d^0$ is the initial depth of the crevice and $C_{\mathrm{Cl}^-}^0$ is given to satisfy the local electroneutrality condition of the solution.

Figure 4.1 shows the solution chemistry obtained by Sharland's steady-state model where the moving boundary is ignored and the pit/crevice bottom position is kept fixed and constant. Figure 4.2 shows the concentration profiles of the ions and the crevice depth calculated by our model at two different times of propagation for the same initial conditions as the ones used in Figure 4.1 from our model.
<figure><img src="image_3.png" /><figcaption>Figure 4.1: Concentration profiles along the cavity length for a crevice with passive walls given by Sharland model [6].</figcaption></figure>

The concentration profiles of ions computed by our diffusion-migration-reaction model [12, p.158] (Figure 4.2) and those by the Sharland et al. steady-state model (Figure 4.1) show a good qualitative agreement. On the other hand, they are not exactly the same quantitatively. This quantitative discrepancy is due to the fact that Sharland et al's model is time-independent. In consequence, it cannot simulate the exact quantitative profiles of the ions in the crevice at each time since the crevice depth evolves with time. Crevice corrosion is always a non steady state process.

<figure><img src="image_4.png" /><figcaption>Figure 4.2: Evolution of the solution chemistry in the crevice given by the developed model with moving interface for the initial data : $\phi_m = 0.04 \text{ V}_{|\text{NHE}}$, $C_{\text{Fe}2+}^0 = 10^{-6} \text{ mol/L}$, $C_{\text{NaCl}}^0 = 10^{-3} \text{ mol/L}$ and $x_d(0) = 2 \text{ mm}$.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 4.1: Concentration profiles along the cavity length for a crevice with passive walls given by Sharland model [6].</figcaption></figure>

The concentration profiles of ions computed by our diffusion-migration-reaction model [12, p.158] (Figure 4.2) and those by the Sharland et al. steady-state model (Figure 4.1) show a good qualitative agreement. On the other hand, they are not exactly the same quantitatively. This quantitative discrepancy is due to the fact that Sharland et al's model is time-independent. In consequence, it cannot simulate the exact quantitative profiles of the ions in the crevice at each time since the crevice depth evolves with time. Crevice corrosion is always a non steady state process.

<figure><img src="image_4.png" /><figcaption>Figure 4.2: Evolution of the solution chemistry in the crevice given by the developed model with moving interface for the initial data : $\phi_m = 0.04 \text{ V}_{\text{NHE}}$, $C_{\text{Fe}2+}^0 = 10^{-6} \text{ mol/L}$, $C_{\text{NaCl}}^0 = 10^{-3} \text{ mol/L}$ and $x_d(0) = 2 \text{ mm}$.</figcaption></figure>
## 4.2. Computation of the minimum value $\phi_m$ for propagation simulation

The electric potential $\phi_m$ applied on the pure iron surface has to be chosen such that the metal initially dissolves ($f_{Fe^{2+}} > 0$) and does not re-deposit (reduction phase corresponding to $f_{Fe^{2+}} < 0$). The precipitation of salt films, often observed in the bottom of the pit or crevice is not considered here. The minimum value of $\phi_m$ to be imposed depends on the initial data. Let $\phi_m^{min}$ be such a value. Then, $\phi_m^{min}$ satisfies the Nernst equation which is given in its dimensional form by:

$$ \varepsilon^{eq} ( := \phi_m^{min} - \phi^* ) = E_{Fe^{2+}/Fe}^0 + \frac{RT}{z_{Fe^{2+}}F} \ln \left( \frac{C_{Fe^{2+}}}{C_{ref}} \right), \quad C_{ref} = 1 \text{ mol/L} \quad (4.2) $$

where $\varepsilon^{eq}$ is the electrode potential of iron in contact with the pit/crevice solution at equilibrium, $E_{Fe^{2+}/Fe}^0$ denotes the standard potential for $Fe^{2+} \leftrightarrow Fe$, which is known to be $E_{Fe^{2+}/Fe}^0 = -0.44 \text{ V}$ relative to the standard hydrogen electrode at room temperature. The value $\phi^*$ is usually called $\phi^{eq}$ in the literature. It denotes the potential at equilibrium for which the Butler-Volmer flux (2.7) vanishes. The curve of $\phi_m^{min} - \phi^*$ as a function of the initial concentration of the iron cations at the bottom of the pit/crevice $C_{Fe^{2+}}^0(x_d^0)$ is plotted in Figure 4.3 (using equation (4.2)).

<figure><img src="image_5.png" /><figcaption>Figure 4.3: $\phi_m^{min} - \phi^*$ as a function of the initial concentration of iron cations at the pit/crevice bottom, given by (4.2).</figcaption></figure>

The region above the red curve corresponds to the range of $\phi_m$ where dissolution occurs ($\phi_m > \phi_m^{min} - \phi^*$ and $f_{Fe^{2+}} > 0$) and the region below corresponds to the reduction field ($\phi_m < \phi_m^{min} - \phi^*$). To be close to equilibrium, the $\phi_m$ value needs to be close to the red curve ($\phi_m \approx \phi_m^{min} - \phi^*$ and $f_{Fe^{2+}} < 0$).

## 4.3. Parametric study: crevice corrosion

We begin with studying crevice corrosion. The pitting corrosion will be discussed in subsection 4.4. The model set up makes it possible to describe the behavior of the system in presence of different variations of the physico-chemical parameters such as:

* the influence of the initial chloride concentration $C_{Cl-}^0$ in the pit/crevice solution.
## 4.2. Computation of the minimum value $\phi_m$ for propagation simulation

The electric potential $\phi_m$ applied on the pure iron surface has to be chosen such that the metal initially dissolves ($f_{Fe^{2+}} > 0$) and does not re-deposit (reduction phase corresponding to $f_{Fe^{2+}} < 0$). The precipitation of salt films, often observed in the bottom of the pit or crevice is not considered here. The minimum value of $\phi_m$ to be imposed depends on the initial data. Let $\phi_m^{min}$ be such a value. Then, $\phi_m^{min}$ satisfies the Nernst equation which is given in its dimensional form by:

$$ \varepsilon^{eq} ( := \phi_m^{min} - \phi^* ) = E_{Fe^{2+}/Fe}^0 + \frac{RT}{z_{Fe^{2+}}F} \ln \left( \frac{C_{Fe^{2+}}}{C_{ref}} \right), \quad C_{ref} = 1 \text{ mol/L} \quad (4.2) $$

where $\varepsilon^{eq}$ is the electrode potential of iron in contact with the pit/crevice solution at equilibrium, $E_{Fe^{2+}/Fe}^0$ denotes the standard potential for $Fe^{2+} \leftrightarrow Fe$, which is known to be $E_{Fe^{2+}/Fe}^0 = -0.44 \text{ V}$ relative to the standard hydrogen electrode at room temperature. The value $\phi^*$ is usually called $\phi^{eq}$ in the literature. It denotes the potential at equilibrium for which the Butler-Volmer flux (2.7) vanishes. The curve of $\phi_m^{min} - \phi^*$ as a function of the initial concentration of the iron cations at the bottom of the pit/crevice $C_{Fe^{2+}}^0(x_d^0)$ is plotted in Figure 4.3 (using equation (4.2)).

<figure><img src="image_5.png" /><figcaption>Figure 4.3: $\phi_m^{min} - \phi^*$ as a function of the initial concentration of iron cations at the pit/crevice bottom, given by (4.2).</figcaption></figure>

The region above the red curve corresponds to the range of $\phi_m$ where dissolution occurs ($\phi_m > \phi_m^{min} - \phi^*$ and $f_{Fe^{2+}} > 0$) and the region below corresponds to the reduction field ($\phi_m < \phi_m^{min} - \phi^*$). To be close to equilibrium, the $\phi_m$ value needs to be close to the red curve ($\phi_m \approx \phi_m^{min} - \phi^*$ and $f_{Fe^{2+}} < 0$).

## 4.3. Parametric study: crevice corrosion

We begin with studying crevice corrosion. The pitting corrosion will be discussed in subsection 4.4. The model set up makes it possible to describe the behavior of the system in presence of different variations of the physico-chemical parameters such as:

* the influence of the initial chloride concentration $C_{Cl-}^0$ in the pit/crevice solution.
* the impact of the potential applied to the surface of pure iron $\phi_m$.

The objective of this subsection is to evaluate the influence of these different parameters on the system without taking into account any passivation criteria. In other words, we do not consider any phenomenon of repassivation if the criterion of crevice stability i.e. the propagation of the cervice (corresponding to a positive velocity displacement of the moving boundary) is not respected. In this context, independently of the value of the applied metal potential $\phi_m$, it is assumed that the crevice will always propagate even for low metallic cation concentration at the bottom of the crevice.

In this subsection, numerical simulations of solutions of system (3.1),(2.22),(2.23) (see Section 3 and Appendix A, Appendix B) for a one-dimensional crevice growth which initially contains $10^{-6}$ mol/L of $Fe^{2+}$, are presented. The computations are done for the following choice of reference parameters

$$D_0 = 1 \cdot 10^{-9} \text{ m}^2/\text{s}, \quad L = 10^3 \text{ µm} = 1 \text{ mm}, \quad C_0 = 10^{-3} \text{ mol/L}, \qquad (4.3)$$

and for the following physical and chemical parameters [23]:

$$D_{Fe^{2+}} = 7.2 \cdot 10^{-10} \text{ m}^2 \cdot \text{s}^{-1}, \quad D_{Cl^{-}} = 2 \cdot 10^{-9} \text{ m}^2 \cdot \text{s}^{-1}, \quad D_{Na^{+}} = 1.3 \cdot 10^{-9} \text{ m}^2 \cdot \text{s}^{-1}, \qquad (4.4)$$

$$k_a = 89.1 \text{ mol.m}^{-2}\text{s}^{-1}, \quad k_c = 1.2 \cdot 10^{-13} \text{ mol.m}^{-2}\text{s}^{-1} \quad (\text{Appendix B of [12, p.195]}) \qquad (4.5)$$

with an initial crevice depth $x_d^0 := x_d(0) = 1$ mm.

### 4.3.1. Influence of the initial chloride concentration $C_{Cl^{-}}^0$.

In this paragraph, simulations are given for a fixed value of $\phi_m = -0.2 \text{ V}_{\text{NHE}}$. Figure 4.4 shows the crevice propagation kinetics for several initial chloride concentrations between 0.001 mol/L and 1 mol/L. It shows that the speed of the crevice propagation increases with the initial chloride concentration. After 500 hours of propagation, the corrosion speed for $C_{NaCl}^0 = 1$ mol/L is roughly 12 times faster than for $C_{NaCl}^0 = 0.001$ mol/L.

<figure><img src="image_6.png" /><figcaption>Figure 4.4: Influence of the initial concentration of chloride ions on the time evolution of the crevice depth.</figcaption></figure>

Figure 4.5 shows the value of the final crevice depth after 500 hours of propagation for several initial chloride concentrations between 0.001 mol/L and 1 mol/L (given by the symbol ♦).
* the impact of the potential applied to the surface of pure iron $\phi_m$.

The objective of this subsection is to evaluate the influence of these different parameters on the system without taking into account any passivation criteria. In other words, we do not consider any phenomenon of repassivation if the criterion of crevice stability i.e. the propagation of the cervice (corresponding to a positive velocity displacement of the moving boundary) is not respected. In this context, independently of the value of the applied metal potential $\phi_m$, it is assumed that the crevice will always propagate even for low metallic cation concentration at the bottom of the crevice.

In this subsection, numerical simulations of solutions of system (3.1),(2.22),(2.23) (see Section 3 and Appendix A, Appendix B) for a one-dimensional crevice growth which initially contains $10^{-6}$ mol/L of $Fe^{2+}$, are presented. The computations are done for the following choice of reference parameters

$$D_0 = 1 \cdot 10^{-9} \text{ m}^2/\text{s}, \quad L = 10^3 \text{ µm} = 1 \text{ mm}, \quad C_0 = 10^{-3} \text{ mol/L}, \qquad (4.3)$$

and for the following physical and chemical parameters [23]:

$$D_{Fe^{2+}} = 7.2 \cdot 10^{-10} \text{ m}^2 \cdot \text{s}^{-1}, \quad D_{Cl^{-}} = 2 \cdot 10^{-9} \text{ m}^2 \cdot \text{s}^{-1}, \quad D_{Na^{+}} = 1.3 \cdot 10^{-9} \text{ m}^2 \cdot \text{s}^{-1}, \qquad (4.4)$$

$$k_a = 89.1 \text{ mol.m}^{-2}\text{s}^{-1}, \quad k_c = 1.2 \cdot 10^{-13} \text{ mol.m}^{-2}\text{s}^{-1} \quad (\text{Appendix B of [12, p.195]}) \qquad (4.5)$$

with an initial crevice depth $x_d^0 := x_d(0) = 1$ mm.

### 4.3.1. Influence of the initial chloride concentration $C_{Cl^-}$.

In this paragraph, simulations are given for a fixed value of $\phi_m = -0.2 \text{ V}_{\text{NHE}}$. Figure 4.4 shows the crevice propagation kinetics for several initial chloride concentrations between 0.001 mol/L and 1 mol/L. It shows that the speed of the crevice propagation increases with the initial chloride concentration. After 500 hours of propagation, the corrosion speed for $C_{NaCl}^0 = 1$ mol/L is roughly 12 times faster than for $C_{NaCl}^0 = 0.001$ mol/L.

<figure><img src="image_6.png" /><figcaption>Figure 4.4: Influence of the initial concentration of chloride ions on the time evolution of the crevice depth.</figcaption></figure>

Figure 4.5 shows the value of the final crevice depth after 500 hours of propagation for several initial chloride concentrations between 0.001 mol/L and 1 mol/L (given by the symbol ♦).
<figure><img src="image_7.png" /><figcaption>Figure 4.5: Evolution of the crevice depth after 500 hours of the crevice propagation for different values of initial chloride concentrations $C_{\text{NaCl}}^0$ (with $x_d^0 = 1$ mm, $\phi_m = -0.2 \text{ V(NHE)}$, $C_{\text{Fe}^{2+}}^0 = 10^{-6} \text{ mol/L}$, $t = 500$ hours).</figcaption></figure>

Since the Butler-Volmer formula (2.7) describes the dissolution rate at the crevice bottom, it follows that this flux incorporates the main factors that influence the corrosion speed. This flux is defined as a function of the iron ions concentration and the potential drop at the crevice bottom. So, in order to physically interpret the results of Figure 4.4, the evolution of the metallic cations concentration $C_{\text{Fe}^{2+}}$ at the crevice bottom after ≈ 500 hours of crevice propagation, is plotted in Figure 4.6 for different values of initial chloride concentrations $C_{\text{NaCl}}^0$.

<figure><img src="image_8.png" /><figcaption>Figure 4.6: Influence of the initial concentration of NaCl on the time evolution of the metallic cation concentration $C_{\text{Fe}^{2+}}$ at the crevice bottom.</figcaption></figure>

Figure 4.6 shows that for $\phi_m = -0.2 \text{ V}|_{\text{NHE}}$, the concentration of $\text{Fe}^{2+}$ at the bottom of the crevice increases until reaching a quasi-stationary value ≈ 82.2 mol/L. This high concentration
<figure><img src="image_7.png" /><figcaption>Figure 4.5: Evolution of the crevice depth after 500 hours of the crevice propagation for different values of initial chloride concentrations $C_{\text{NaCl}}^0$ (with $x_d^0 = 1$ mm, $\phi_m = -0.2 \text{ V(NHE)}$, $C_{\text{Fe}^{2+}}^0 = 10^{-6} \text{ mol/L}, t = 500$ hours).</figcaption></figure>

Since the Butler-Volmer formula (2.7) describes the dissolution rate at the crevice bottom, it follows that this flux incorporates the main factors that influence the corrosion speed. This flux is defined as a function of the iron ions concentration and the potential drop at the crevice bottom. So, in order to physically interpret the results of Figure 4.4, the evolution of the metallic cations concentration $C_{\text{Fe}^{2+}}$ at the crevice bottom after $\approx 500$ hours of crevice propagation, is plotted in Figure 4.6 for different values of initial chloride concentrations $C_{\text{NaCl}}^0$.

<figure><img src="image_8.png" /><figcaption>Figure 4.6: Influence of the initial concentration of NaCl on the time evolution of the metallic cation concentration $C_{\text{Fe}^{2+}}$ at the crevice bottom.</figcaption></figure>

Figure 4.6 shows that for $\phi_m = -0.2 \text{ V}|_{\text{NHE}}$, the concentration of $\text{Fe}^{2+}$ at the bottom of the crevice increases until reaching a quasi-stationary value $\approx 82.2 \text{ mol/L}$. This high concentration
value was also reached for lower initial chloride concentrations $C_{\text{NaCl}}^0$ ($< 1 \text{ mol/L}$) but after much longer times. The transient time required to reach this near-equilibrium state decreases when $C_{\text{NaCl}}^0$ increases. It is obvious that the situation where very high $\text{Fe}^{2+}$ concentration are reached does not reflect reality because these iron ions precipitate into chloride salt film if the metallic cations concentration $C_{\text{Fe}^{2+}}$ exceeds the saturation value of 5 mol/L at the crevice bottom.

Figure 4.7 describes the evolution of the potential drop $\Delta\phi = \phi(x_d(t), t) - \phi(0, t) = \phi(x_d(t), t)$, which coincides, in this case, with the potential at the bottom of the crevice $\phi(x_d(t), t)$. It shows that the potential drop converges to a certain constant potential value that depends on $C_{\text{NaCl}}^0$. Physically, this can be described as a convergence to a state of equilibrium after a large time of crevice propagation : $C_{\text{Fe}^{2+}}(x_d(t), t) \approx$ constant and $\phi(x_d(t), t) \approx$ constant as $t \to \infty$. Incidentally, as we see in Figure 4.7, the lower the NaCl concentration, the higher the potential drop. This is quite natural since lower ionic concentrations imply higher electric resistance of the solution.

<figure><img src="image_9.png" /><figcaption>Figure 4.7: Influence of the initial concentration of NaCl on the time evolution of the potential drop in the crevice solution with time.</figcaption></figure>

Interestingly, the potential drop for the higher NaCl concentrations increases and then decreases. This behavior could be induced by the strong increase of the iron concentration observed in Figure 4.6 which increases the solution conductivity faster than the pit depth growth.

### 4.3.2. Influence of the iron potential $\phi_m$

The simulations are given for a fixed value of the initial chloride concentration $C_{\text{NaCl}}^0 = 10^{-3} \text{ mol/L}$. Figure 4.8 shows the graphs of the time evolution of crevice depths for different values of $\phi_m$ and Figure 4.9, the final crevice depth after 500 hours of propagation for different values of the metal potential between $-0.35\text{V}|_{\text{NHE}}$ and $-0.2\text{V}|_{\text{NHE}}$.

Both graphs illustrate the very high sensitivity of the propagation speed to the metal potential value $\phi_m$. Figure 4.8 shows that the evolution of the corrosion speed is constant for the lower values of $\phi_m$. Moreover, it is clear that for the higher values of $\phi_m$ (for example, $\phi_m = -0.2 \text{ V}$), the crevice cannot keep a constant speed of growth.
value was also reached for lower initial chloride concentrations $C_{\text{NaCl}}^0$ ($<1$ mol/L) but after much longer times. The transient time required to reach this near-equilibrium state decreases when $C_{\text{NaCl}}^0$ increases. It is obvious that the situation where very high $\text{Fe}^{2+}$ concentration are reached does not reflect reality because these iron ions precipitate into chloride salt film if the metallic cations concentration $C_{\text{Fe}^{2+}}$ exceeds the saturation value of 5 mol/L at the crevice bottom.

Figure 4.7 describes the evolution of the potential drop $\Delta\phi = \phi(x_d(t), t) - \phi(0, t) = \phi(x_d(t), t)$, which coincides, in this case, with the potential at the bottom of the crevice $\phi(x_d(t), t)$. It shows that the potential drop converges to a certain constant potential value that depends on $C_{\text{NaCl}}^0$. Physically, this can be described as a convergence to a state of equilibrium after a large time of crevice propagation : $C_{\text{Fe}^{2+}}(x_d(t), t) \approx$ constant and $\phi(x_d(t), t) \approx$ constant as $t \to \infty$. Incidentally, as we see in Figure 4.7, the lower the NaCl concentration, the higher the potential drop. This is quite natural since lower ionic concentrations imply higher electric resistance of the solution.

<figure><img src="image_9.png" /><figcaption>Figure 4.7: Influence of the initial concentration of NaCl on the time evolution of the potential drop in the crevice solution with time.</figcaption></figure>

Interestingly, the potential drop for the higher NaCl concentrations increases and then decreases. This behavior could be induced by the strong increase of the iron concentration observed in Figure 4.6 which increases the solution conductivity faster than the pit depth growth.

### 4.3.2. Influence of the iron potential $\phi_m$

The simulations are given for a fixed value of the initial chloride concentration $C_{\text{NaCl}}^0 = 10^{-3}$ mol/L. Figure 4.8 shows the graphs of the time evolution of crevice depths for different values of $\phi_m$ and Figure 4.9, the final crevice depth after 500 hours of propagation for different values of the metal potential between $-0.35\text{V}|_{\text{NHE}}$ and $-0.2\text{V}|_{\text{NHE}}$.

Both graphs illustrate the very high sensitivity of the propagation speed to the metal potential value $\phi_m$. Figure 4.8 shows that the evolution of the corrosion speed is constant for the lower values of $\phi_m$. Moreover, it is clear that for the higher values of $\phi_m$ (for example, $\phi_m = -0.2$ V), the crevice cannot keep a constant speed of growth.
<figure><img src="image_10.png" /><figcaption>Figure 4.8: Influence of the potential applied to the surface of the pure iron $\phi_m$ on the evolution of the crevice depth in time.</figcaption></figure>

Figure 4.9 shows the value of the final crevice depth after 500 hours of propagation for different values of the metal potential $\phi_m$ between $-0.35 \text{ V}|_{\text{NHE}}$ and $-0.2 \text{ V}|_{\text{NHE}}$.

<figure><img src="image_11.png" /><figcaption>Figure 4.9: Estimation of the depth evolution law after 500 hours of the crevice propagation as a function of the metal potential $\phi_m$.</figcaption></figure>

Figure 4.10 and Figure 4.11 show the evolution of the iron concentration at the crevice bottom and the evolution of the potential drop with time, respectively. It is observed that both parameters, after an initial very fast transient stage, are quite constant with time for the lower potentials but increases with time for the higher potentials. Since the corrosion rate driven by the Butler-
<figure><img src="image_10.png" /><figcaption>Figure 4.8: Influence of the potential applied to the surface of the pure iron $\phi_m$ on the evolution of the crevice depth in time.</figcaption></figure>

Figure 4.9 shows the value of the final crevice depth after 500 hours of propagation for different values of the metal potential $\phi_m$ between $-0.35 \text{ V}|_{\text{NHE}}$ and $-0.2 \text{ V}|_{\text{NHE}}$.

<figure><img src="image_11.png" /><figcaption>Figure 4.9: Estimation of the depth evolution law after 500 hours of the crevice propagation as a function of the metal potential $\phi_m$.</figcaption></figure>

Figure 4.10 and Figure 4.11 show the evolution of the iron concentration at the crevice bottom and the evolution of the potential drop with time, respectively. It is observed that both parameters, after an initial very fast transient stage, are quite constant with time for the lower potentials but increases with time for the higher potentials. Since the corrosion rate driven by the Butler-
Volmer formula is dependent on both parameters, their evolution explains the constant or increasing corrosion rate observed in Figure 4.8.

<figure><img src="image_12.png" /><figcaption>Figure 4.10: Influence of the metal potential $\phi_m$ on the time evolution of the metallic cation concentration $C_{Fe^{2+}}$ at the crevice bottom.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 4.11: Influence of the metal potential $\phi_m$ on the time evolution of the potential drop in the crevice solution with time.</figcaption></figure>
Volmer formula is dependent on both parameters, their evolution explains the constant or increasing corrosion rate observed in Figure 4.8.

<figure><img src="image_12.png" /><figcaption>Figure 4.10: Influence of the metal potential $\phi_m$ on the time evolution of the metallic cation concentration $C_{Fe^{2+}}$ at the crevice bottom.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 4.11: Influence of the metal potential $\phi_m$ on the time evolution of the potential drop in the crevice solution with time.</figcaption></figure>
### 4.3.3. Proportion of diffusion and migration during crevice propagation

It is important to determine the proportion of transport by migration and by diffusion in order to evaluate the validity of neglecting the migration term in the governing transport equation. In the literature, some authors consider only transport by diffusion in their calculations to make the equations more easily solvable, but this could lead to significant errors in longtime prediction if the effect of transport by migration is not so small. Therefore, quantifying this error is of vital importance. The mass-conservation equation (in its dimensionless form) that governs the concentrations of aqueous species in the dilute solution filling the pit (or crevice) is given by:

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right), \quad t > 0, \quad 0 < x < x_d(t), \quad j = 1, \dots, M \quad (4.6) $$

The first term on the right hand side of (4.6) describes the rate of transport of ions by diffusion under concentration gradients. The second term represents the migration of charged species under electrostatic potential gradients. Therefore, the ionic flux carried by the j-th ion species, is given in its dimensionless form by:

$$ J_j := J_j^{\text{diff}} + J_j^{\text{mig}} = -D_j \frac{\partial C_j}{\partial x} - D_j z_j C_j \frac{\partial \phi}{\partial x}, \quad (4.7) $$

where $J_j^{\text{diff}} = -D_j \frac{\partial C_j}{\partial x}$ is the diffusion flux which is associated with Fick's law and $J_j^{\text{mig}} = -D_j z_j C_j \frac{\partial \phi}{\partial x}$ is the migration flux caused by the electric field $E = -\frac{\partial \phi}{\partial x}$. The proportion of transport by migration of the species j versus total transport is computed as follows:

$$ \% \text{ migration\_transport\_j} = 100 \cdot \frac{|J_j^{\text{mig}}|}{|J_j^{\text{diff}}| + |J_j^{\text{mig}}|} \quad (4.8) $$

Two cases of localized corrosion were evaluated: crevice corrosion and pit corrosion. For crevice corrosion, an initial crevice depth $x_d^0 = 100 \mu\text{m}$ and an initial iron concentration $C_{\text{Fe}^{2+}}^0 = 10^{-6} \text{mol/L}$ were assumed according to usual crevice configuration. The computations showed that the transport of $\text{Cl}^-$ and $\text{Na}^+$ ions is about 50% by diffusion and 50% by migration along the pit solution, regardless of the metal potential value and the initial concentration of NaCl in the solution. However, the results for the $\text{Fe}^{2+}$ transport were quite different. Figures 4.12 and 4.13 show the profile of the $\text{Fe}^{2+}$ concentration and that of the potential at each point of the crevice solution, for different values of NaCl concentration and the metal potential $\phi_m$ when the crevice depth reaches 200 $\mu\text{m}$. These figures show that the $\text{Fe}^{2+}$ concentration and the potential drop both increase as $\phi_m$ increases. And, for a fixed value of $\phi_m$, lower NaCl concentration leads to lower $\text{Fe}^{2+}$ concentration (Figure 4.12) and higher potential drop (Figure 4.13). This higher potential drop is no surprise, since lower NaCl and $\text{Fe}^{2+}$ concentrations imply higher electric resistance of the crevice solution.

Figures 4.14 and 4.15 show the proportion of transport by migration for $\text{Fe}^{2+}$ species for the same set of values of NaCl concentration and $\phi_m$ as in Figures 4.12 and 4.13, when the crevice depth reaches 200 $\mu\text{m}$ (Figure 4.14) and 1 mm (Figure 4.15). Regardless of the initial NaCl concentration varying between 0.001 mol/L and 1 mol/L, the proportion of $\text{Fe}^{2+}$ transport by migration increases with the metal potential value $\phi_m$. Moreover, for a fixed value of $\phi_m$, the proportion becomes more significant for lower initial concentrations of NaCl. Figures 4.14 and 4.15 demonstrate that the proportion of migration of $\text{Fe}^{2+}$ cations in the crevice solution becomes substantially higher when
### 4.3.3. Proportion of diffusion and migration during crevice propagation

It is important to determine the proportion of transport by migration and by diffusion in order to evaluate the validity of neglecting the migration term in the governing transport equation. In the literature, some authors consider only transport by diffusion in their calculations to make the equations more easily solvable, but this could lead to significant errors in longtime prediction if the effect of transport by migration is not so small. Therefore, quantifying this error is of vital importance. The mass-conservation equation (in its dimensionless form) that governs the concentrations of aqueous species in the dilute solution filling the pit (or crevice) is given by:

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right), \quad t > 0, \quad 0 < x < x_d(t), \quad j = 1, \dots, M \quad (4.6) $$

The first term on the right hand side of (4.6) describes the rate of transport of ions by diffusion under concentration gradients. The second term represents the migration of charged species under electrostatic potential gradients. Therefore, the ionic flux carried by the j-th ion species, is given in its dimensionless form by:

$$ J_j := J_j^{\text{diff}} + J_j^{\text{mig}} = -D_j \frac{\partial C_j}{\partial x} - D_j z_j C_j \frac{\partial \phi}{\partial x}, \quad (4.7) $$

where $J_j^{\text{diff}} = -D_j \frac{\partial C_j}{\partial x}$ is the diffusion flux which is associated with Fick's law and $J_j^{\text{mig}} = -D_j z_j C_j \frac{\partial \phi}{\partial x}$ is the migration flux caused by the electric field $E = -\frac{\partial \phi}{\partial x}$. The proportion of transport by migration of the species j versus total transport is computed as follows:

$$ \% \text{ migration\_transport\_j} = 100 \cdot \frac{|J_j^{\text{mig}}|}{|J_j^{\text{diff}}| + |J_j^{\text{mig}}|} \quad (4.8) $$

Two cases of localized corrosion were evaluated: crevice corrosion and pit corrosion. For crevice corrosion, an initial crevice depth $x_d^0 = 100 \mu\text{m}$ and an initial iron concentration $C_{\text{Fe}^{2+}}^0 = 10^{-6} \text{mol/L}$ were assumed according to usual crevice configuration. The computations showed that the transport of $\text{Cl}^-$ and $\text{Na}^+$ ions is about 50% by diffusion and 50% by migration along the pit solution, regardless of the metal potential value and the initial concentration of NaCl in the solution. However, the results for the $\text{Fe}^{2+}$ transport were quite different. Figures 4.12 and 4.13 show the profile of the $\text{Fe}^{2+}$ concentration and that of the potential at each point of the crevice solution, for different values of NaCl concentration and the metal potential $\phi_m$ when the crevice depth reaches 200 $\mu\text{m}$. These figures show that the $\text{Fe}^{2+}$ concentration and the potential drop both increase as $\phi_m$ increases. And, for a fixed value of $\phi_m$, lower NaCl concentration leads to lower $\text{Fe}^{2+}$ concentration (Figure 4.12) and higher potential drop (Figure 4.13). This higher potential drop is no surprise, since lower NaCl and $\text{Fe}^{2+}$ concentrations imply higher electric resistance of the crevice solution.

Figures 4.14 and 4.15 show the proportion of transport by migration for $\text{Fe}^{2+}$ species for the same set of values of NaCl concentration and $\phi_m$ as in Figures 4.12 and 4.13, when the crevice depth reaches 200 $\mu\text{m}$ (Figure 4.14) and 1 mm (Figure 4.15). Regardless of the initial NaCl concentration varying between 0.001 mol/L and 1 mol/L, the proportion of $\text{Fe}^{2+}$ transport by migration increases with the metal potential value $\phi_m$. Moreover, for a fixed value of $\phi_m$, the proportion becomes more significant for lower initial concentrations of NaCl. Figures 4.14 and 4.15 demonstrate that the proportion of migration of $\text{Fe}^{2+}$ cations in the crevice solution becomes substantially higher when
the NaCl concentrations are smaller. As shown in Figures 4.12 and 4.13, a lower NaCl concentration leads to a lower $Fe^{2+}$ concentration and a larger potential drop. Apparently, such an environment promotes higher proportion of migration of $Fe^{2+}$ cations in the crevice solution.

<figure><img src="image_14.png" /><figcaption>Figure 4.12: Profile of $Fe^{2+}$ concentration in crevice solution as function of metal potential and NaCl concentration.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 4.13: Potential drop in crevice solution as function of metal potential and NaCl concentration.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 4.14: Proportion of migration in transport (%) in crevice solution as function of metal potential and NaCl concentration for a 200µm deep crevice.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 4.15: Proportion of migration in transport (%) in crevice solution as function of metal potential and NaCl concentration for a 1 mm deep crevice.</figcaption></figure>

Transport by migration is considered to be negligible only in the case of a low metal potential, -0.4 V/NHE with high initial concentration of NaCl, 1mol/L. Such cases could be encountered in pit configuration but needs to be checked. It is done below.

### 4.4. Parametric study for pit propagation

#### 4.4.1. Evaluation of the critical conditions for pit propagation

Pit configuration differs from crevice corrosion in that pit propagation often starts at a very shallow depth, typically a few to hundreds of nanometers. Additionally, a salt film with a high
the NaCl concentrations are smaller. As shown in Figures 4.12 and 4.13, a lower NaCl concentration leads to a lower $Fe^{2+}$ concentration and a larger potential drop. Apparently, such an environment promotes higher proportion of migration of $Fe^{2+}$ cations in the crevice solution.

<figure><img src="image_14.png" /><figcaption>Figure 4.12: Profile of $Fe^{2+}$ concentration in crevice solution as function of metal potential and NaCl concentration.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 4.13: Potential drop in crevice solution as function of metal potential and NaCl concentration.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 4.14: Proportion of migration in transport (%) in crevice solution as function of metal potential and NaCl concentration for a 200µm deep crevice.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 4.15: Proportion of migration in transport (%) in crevice solution as function of metal potential and NaCl concentration for a 1 mm deep crevice.</figcaption></figure>

Transport by migration is considered to be negligible only in the case of a low metal potential, -0.4 V/NHE with high initial concentration of NaCl, 1mol/L. Such cases could be encountered in pit configuration but needs to be checked. It is done below.

### 4.4. Parametric study for pit propagation

#### 4.4.1. Evaluation of the critical conditions for pit propagation

Pit configuration differs from crevice corrosion in that pit propagation often starts at a very shallow depth, typically a few to hundreds of nanometers. Additionally, a salt film with a high
concentration of chloride ions is present on the metallic surface in contact with the pit solution [24]. It is commonly proposed in the literature that a critical concentration of metal cation at the bottom of the pit, denoted by $C_{crit}$, is necessary to maintain pit propagation. The maintenance of this critical condition at the pit bottom depends heavily on the metal potential and the geometry of the pit, such as its depth. These dependencies have been computed using our model. In the simulations that follow, the critical concentration of metal cation is assumed to be about 60% of the saturation of the metallic cation, as proposed in [24]. Therefore, to ensure the stability of the pit, the metallic cation concentration at the pit bottom $C_{Fe^{2+}}(x_d(t), t)$ should be greater than $C_{crit} = 60\% C_{sat}$ at any time. Laycock et al. [24] assumed that for the 300 series alloy model, the saturated concentration $C_{sat}$ is equal to 5 mol/L. Srinivasan et al. [25] used a saturated concentration of the 316L cation equal to 5.02 mol/L. Thus, in the case of pure iron, $C_{sat}$ is assumed to be 5 mol/L, and $C_{crit}$ is equal to 3 mol/L. A classical pit propagation scenario begins with a pit nuclei that is covered by a lacy cover, which helps to maintain a concentrated local chemistry [24]. In fact, many pits are entirely dependent on the existence of these covers for their continued propagation. However, a long-term propagating pit has the particularity of being able to propagate without any lacy cover. Therefore, the objective of the following simulations is to identify the required conditions on the metal potential $\phi_m$, to maintain the stability of the pit without the protection of these lacy covers. The critical potential for maintaining pit propagation, once the lacy covers have collapsed, strongly depends on the reached pit depth [24]. For these simulations, the following system is considered: at the pit entrance ($x = 0$) and at time $t = 0$, the concentration of each aqueous species is equal to its value in the bulk solution. At this time $t = 0$, a salt film is considered to have formed on the bottom of the pit, implying that the concentration of iron cations has reached its saturation value, $C_{Fe^{2+}}(x_d(0), 0) = 5$ mol/L. The following choice of initial data is used for simulations:

At the pit entrance ($x = 0$ and at time $t = 0$):

$$C_{Fe^{2+}}(0, 0) = 10^{-6} \text{ mol/L} \quad \text{and} \quad C_{Na^{+}}(0, 0) = C_{Cl^{-}}(0, 0) = 1 \text{ mol/L.} \qquad (4.9)$$

At the pit bottom ($x = x_d(0)$):

$$\begin{aligned} C_{Fe^{2+}}(x_d(0), 0) &= 5 \text{ mol/L, } C_{Na^{+}}(x_d(0), 0) = 1 \text{ mol/L,} \\ C_{Cl^{-}}(x_d(0), 0) &= 11 \text{ mol/L (by imposed local electroneutrality of the solution).} \end{aligned} \qquad (4.10)$$

In the pit solution ($0 < x < x_d(0)$):

$$\begin{aligned} C_{Fe^{2+}}(x, 0) &= \text{linear profile from } 10^{-6} \text{ mol/L at } x = 0 \text{ to } 5 \text{ mol/L at } x = x_d(0), \\ C_{Cl^{-}}(x, 0) &= \text{linear profile from } 1 \text{ mol/L at } x = 0 \text{ to } 11 \text{ mol/L at } x = x_d(0), \\ C_{Na^{+}}(x, 0) &= 1 \text{ mol/L for all } 0 < x < x_d(0) \end{aligned} \qquad (4.11)$$

The repassivation potential $\phi_m^{stable}$, is a crucial parameter that defines the potential at which the concentration of iron at the pit bottom falls below 3 mol/L, the critical iron concentration for pit stability. In Figure 4.16, the repassivation potential $\phi_m^{stable}$ is shown relative to the NHE reference for several initial pit depths $x_d^0 := x_d(0)$. As observed, the repassivation potential $\phi_m^{stable}$ is a decreasing function of the initial pit depth. This finding is consistent with the experimental repassivation potential shown in Figure 4.17 (see [25, 26]).
concentration of chloride ions is present on the metallic surface in contact with the pit solution [24]. It is commonly proposed in the literature that a critical concentration of metal cation at the bottom of the pit, denoted by $C_{crit}$, is necessary to maintain pit propagation. The maintenance of this critical condition at the pit bottom depends heavily on the metal potential and the geometry of the pit, such as its depth. These dependencies have been computed using our model. In the simulations that follow, the critical concentration of metal cation is assumed to be about 60% of the saturation of the metallic cation, as proposed in [24]. Therefore, to ensure the stability of the pit, the metallic cation concentration at the pit bottom $C_{Fe^{2+}}(x_d(t), t)$ should be greater than $C_{crit} = 60\% C_{sat}$ at any time. Laycock et al. [24] assumed that for the 300 series alloy model, the saturated concentration $C_{sat}$ is equal to 5 mol/L. Srinivasan et al. [25] used a saturated concentration of the 316L cation equal to 5.02 mol/L. Thus, in the case of pure iron, $C_{sat}$ is assumed to be 5 mol/L, and $C_{crit}$ is equal to 3 mol/L. A classical pit propagation scenario begins with a pit nuclei that is covered by a lacy cover, which helps to maintain a concentrated local chemistry [24]. In fact, many pits are entirely dependent on the existence of these covers for their continued propagation. However, a long-term propagating pit has the particularity of being able to propagate without any lacy cover. Therefore, the objective of the following simulations is to identify the required conditions on the metal potential $\phi_m$, to maintain the stability of the pit without the protection of these lacy covers. The critical potential for maintaining pit propagation, once the lacy covers have collapsed, strongly depends on the reached pit depth [24]. For these simulations, the following system is considered: at the pit entrance ($x = 0$) and at time $t = 0$, the concentration of each aqueous species is equal to its value in the bulk solution. At this time $t = 0$, a salt film is considered to have formed on the bottom of the pit, implying that the concentration of iron cations has reached its saturation value, $C_{Fe^{2+}}(x_d(0), 0) = 5$ mol/L. The following choice of initial data is used for simulations:

At the pit entrance ($x = 0$ and at time $t = 0$):

$$C_{Fe^{2+}}(0, 0) = 10^{-6} \text{ mol/L} \quad \text{and} \quad C_{Na^{+}}(0, 0) = C_{Cl^{-}}(0, 0) = 1 \text{ mol/L}. \qquad (4.9)$$

At the pit bottom ($x = x_d(0)$):

$$\begin{aligned} C_{Fe^{2+}}(x_d(0), 0) &= 5 \text{ mol/L}, & C_{Na^{+}}(x_d(0), 0) &= 1 \text{ mol/L}, \\ C_{Cl^{-}}(x_d(0), 0) &= 11 \text{ mol/L} \quad (\text{by imposed local electroneutrality of the solution}). \end{aligned} \qquad (4.10)$$

In the pit solution ($0 < x < x_d(0)$):

$$\begin{aligned} C_{Fe^{2+}}(x, 0) &= \text{linear profile from } 10^{-6} \text{ mol/L at } x = 0 \text{ to } 5 \text{ mol/L at } x = x_d(0), \\ C_{Cl^{-}}(x, 0) &= \text{linear profile from } 1 \text{ mol/L at } x = 0 \text{ to } 11 \text{ mol/L at } x = x_d(0), \\ C_{Na^{+}}(x, 0) &= 1 \text{ mol/L for all } 0 < x < x_d(0) \end{aligned} \qquad (4.11)$$

The repassivation potential $\phi_m^{stable}$, is a crucial parameter that defines the potential at which the concentration of iron at the pit bottom falls below 3 mol/L, the critical iron concentration for pit stability. In Figure 4.16, the repassivation potential $\phi_m^{stable}$ is shown relative to the NHE reference for several initial pit depths $x_d^0 := x_d(0)$. As observed, the repassivation potential $\phi_m^{stable}$ is a decreasing function of the initial pit depth. This finding is consistent with the experimental repassivation potential shown in Figure 4.17 (see [25, 26]).
<figure><img src="image_18.png" /><figcaption>Figure 4.16: Identification of the minimal applied potential to ensure the stability of the pit as a function of its initial depth.</figcaption></figure>

<caption>Figure 4.17: Observed behavior of the repassivation potential as a function of the pit depth from artificial pit experiments on stainless steels in chloride media [25].</caption>

The results in Figure 4.16 suggest that the metal potential needs to increase with pit depth to maintain its propagation. For example, for a 10 nm deep pit, the metal potential must be lower than 92 mV to repassivate, whereas for a 100 nm deep pit, it must be lower than 33 mV. This dependence of the repassivation potential on pit depth could explain the wide range of experimental repassivation potential measurements found in literature, as they are most likely carried out at varying initial pit depths. Overall, understanding the relationship between repassivation potential and pit depth is critical to preventing and mitigating pitting corrosion in metals.

### 4.4.2. Proportion of diffusion and migration during pit propagation

The proportion of transport by migration for Fe²⁺ cations was evaluated in the context of an initial pit embryo depth of $x_d(0) = 10$ nm. To define the pit configuration, concentration profiles of Fe²⁺, Cl⁻ and Na⁺ in the pit embryo were used as previously proposed. Based on Figure 4.3, $\phi_m$ is set at 0.1 V/NHE to be higher than the critical $\phi_m$ for an initial pit depth of 10 nm.

The concentration profiles of Fe²⁺, potential profile and proportion of transport by migration in the pit are depicted in Figure 4.18 to 4.20, for a final pit depth of 100 nm. These figures indicate that the proportion of migration for Fe²⁺ transport is quite high, above 60%, at all positions within the pit solution except at the pit mouth, where it decreases significantly. In a pit that has reached 10 µm, the proportion was even higher, approximately 70% (not shown). These observations demonstrate that in conditions allowing pit propagation, the transport of Fe²⁺ by migration cannot be neglected for stability prediction or propagation kinetics, even if the pit solution is highly concentrated and conductive.
<figure><img src="image_18.png" /><figcaption>Figure 4.16: Identification of the minimal applied potential to ensure the stability of the pit as a function of its initial depth.</figcaption></figure>

<caption>Figure 4.17: Observed behavior of the repassivation potential as a function of the pit depth from artificial pit experiments on stainless steels in chloride media [25].</caption>

The results in Figure 4.16 suggest that the metal potential needs to increase with pit depth to maintain its propagation. For example, for a 10 nm deep pit, the metal potential must be lower than 92 mV to repassivate, whereas for a 100 nm deep pit, it must be lower than 33 mV. This dependence of the repassivation potential on pit depth could explain the wide range of experimental repassivation potential measurements found in literature, as they are most likely carried out at varying initial pit depths. Overall, understanding the relationship between repassivation potential and pit depth is critical to preventing and mitigating pitting corrosion in metals.

### 4.4.2. Proportion of diffusion and migration during pit propagation

The proportion of transport by migration for $Fe^{2+}$ cations was evaluated in the context of an initial pit embryo depth of $x_d(0) = 10$ nm. To define the pit configuration, concentration profiles of $Fe^{2+}$, $Cl^-$ and $Na^+$ in the pit embryo were used as previously proposed. Based on Figure 4.3, $\phi_m$ is set at 0.1 V/NHE to be higher than the critical $\phi_m$ for an initial pit depth of 10 nm.

The concentration profiles of $Fe^{2+}$, potential profile and proportion of transport by migration in the pit are depicted in Figure 4.18 to 4.20, for a final pit depth of 100 nm. These figures indicate that the proportion of migration for $Fe^{2+}$ transport is quite high, above 60%, at all positions within the pit solution except at the pit mouth, where it decreases significantly. In a pit that has reached 10 µm, the proportion was even higher, approximately 70% (not shown). These observations demonstrate that in conditions allowing pit propagation, the transport of $Fe^{2+}$ by migration cannot be neglected for stability prediction or propagation kinetics, even if the pit solution is highly concentrated and conductive.
<figure><img src="image_19.png" /><figcaption>Figure 4.18: Profile of $Fe^{2+}$ concentration in 100 nm deep pit.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 4.19: Potential profile in 100 nm deep pit.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 4.20: Proportion of migration in transport (%) in 100 nm deep pit.</figcaption></figure>

To better confirm this conclusion, the kinetics of pit propagation was simulated with a pure diffusion model and compared to the kinetics obtained with the diffusion-migration model presented in this article (see Figure 4.21). The pure diffusion model is composed, on one hand, by the system (2.13) for the concentrations where the migration terms involving derivatives of the potential $\phi$ have been removed (but where of course the potential is kept in the Butler-Volmer flux $f_{Fe^{2+}}(\phi, C_1)$), and on the other hand, by the Partial Differential Equation System (2.12) for the potential in place of the electroneutrality condition (3.1). As a consequence, it must be emphasized that the pure diffusion model ignores the electroneutrality condition which makes this model physically not relevant. The pure diffusion system is numerically solved by decoupling the computation of the potential in (2.12) with the computation of the concentrations in the pure diffusion Partial
<figure><img src="image_19.png" /><figcaption>Figure 4.18: Profile of Fe²⁺ concentration in 100 nm deep pit.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 4.19: Potential profile in 100 nm deep pit.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 4.20: Proportion of migration in transport (%) in 100 nm deep pit.</figcaption></figure>

To better confirm this conclusion, the kinetics of pit propagation was simulated with a pure diffusion model and compared to the kinetics obtained with the diffusion-migration model presented in this article (see Figure 4.21). The pure diffusion model is composed, on one hand, by the system (2.13) for the concentrations where the migration terms involving derivatives of the potential $\phi$ have been removed (but where of course the potential is kept in the Butler-Volmer flux $f_{Fe^{2+}}(\phi, C_1)$), and on the other hand, by the Partial Differential Equation System (2.12) for the potential in place of the electroneutrality condition (3.1). As a consequence, it must be emphasized that the pure diffusion model ignores the electroneutrality condition which makes this model physically not relevant. The pure diffusion system is numerically solved by decoupling the computation of the potential in (2.12) with the computation of the concentrations in the pure diffusion Partial
Differential Equation System. The concentrations are computed at a current time step $t^{n+1}$ by using an semi-implicit scheme (in time) with an explicit term for the Buttler-Volmer flux $f_{Fe^{2+}}(\phi, C_1)$) at the previous time step $t^n$. Then, we compute at time $t^{n+1}$, the profile of the potential $\phi$ at the pit bottom $x = x_d(t^{n+1})$, from known values of concentrations at time $t^{n+1}$. This is done with the calculation explained in Section 3.2.

The comparison between the two models is done for a pit configuration at high NaCl concentration (1 M) with an initial pit depth of 1 µm at metal potential far from equilibrium (-0.2 V/NHE).

<figure><img src="image_22.png" /><figcaption>Figure 4.21: Evolution in time of the pit depth for the pure diffusion model (in yellow) and the diffusion-migration model (in black) for $\phi_m = -0.2$ V/NHE.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 4.22: Potential profiles in the pit obtained for the pure diffusion model (in yellow) and the diffusion-migration model (in black) for $\phi_m = -0.2$ V/NHE.</figcaption></figure>

It is observed in Figure 4.21 that the pure diffusion model overestimates the propagation rate of the pit by a factor of 2 roughly. Thus, migration slows down the pit growth rate. The higher growth rate of the pit when only diffusion is considered is explained by a much lower potential drop in the pit (see Figure 4.22). This discrepancy between simulations from pure diffusion and diffusion-migration models decreases when the overpotential decreases. At a metal potential close to equilibrium, $\phi_m = -0.44$ V/NHE, the pit growth kinetics is roughly equal (not shown). Indeed, the pit reaches a depth of 100 µm after 1194 h for both models. Thus, the discrepancy of the results obtained from the pure diffusion model and the more realistic diffusion-migration model is strongly dependent on the metal potential. In conclusion, the diffusion-migration model must be used for more reliable results at any metal potential.

## 5. Conclusion

A one-dimensional anodic dissolution model with time-dependent diffusion-migration transport was developed for simulating localized corrosion. The mathematical formulation of this physical phenomenon and its computational formulation using a numerical scheme were described in detail to guide corrosion specialists. From a mathematical standpoint, the problem can be identified as a Stefan problem involving a convection-diffusion system. The model was applied to iron with the Butler-Volmer formula as the dissolution law, and both crevice and pit configurations were simulated. The results for crevice corrosion are qualitatively in good agreement with the predictions from the stationary crevice model made by Sharland, but our time-dependent model provides more
Differential Equation System. The concentrations are computed at a current time step $t^{n+1}$ by using an semi-implicit scheme (in time) with an explicit term for the Buttler-Volmer flux $f_{Fe^{2+}}(\phi, C_1)$) at the previous time step $t^n$. Then, we compute at time $t^{n+1}$, the profile of the potential $\phi$ at the pit bottom $x = x_d(t^{n+1})$, from known values of concentrations at time $t^{n+1}$. This is done with the calculation explained in Section 3.2.

The comparison between the two models is done for a pit configuration at high NaCl concentration (1 M) with an initial pit depth of 1 µm at metal potential far from equilibrium (-0.2 V/NHE).

<figure><img src="image_22.png" /><figcaption>Figure 4.21: Evolution in time of the pit depth for the pure diffusion model (in yellow) and the diffusion-migration model (in black) for $\phi_m = -0.2$ V/NHE.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 4.22: Potential profiles in the pit obtained for the pure diffusion model (in yellow) and the diffusion-migration model (in black) for $\phi_m = -0.2$ V/NHE.</figcaption></figure>

It is observed in Figure 4.21 that the pure diffusion model overestimates the propagation rate of the pit by a factor of 2 roughly. Thus, migration slows down the pit growth rate. The higher growth rate of the pit when only diffusion is considered is explained by a much lower potential drop in the pit (see Figure 4.22). This discrepancy between simulations from pure diffusion and diffusion-migration models decreases when the overpotential decreases. At a metal potential close to equilibrium, $\phi_m = -0.44$ V/NHE, the pit growth kinetics is roughly equal (not shown). Indeed, the pit reaches a depth of 100 µm after 1194 h for both models. Thus, the discrepancy of the results obtained from the pure diffusion model and the more realistic diffusion-migration model is strongly dependent on the metal potential. In conclusion, the diffusion-migration model must be used for more reliable results at any metal potential.

## 5. Conclusion

A one-dimensional anodic dissolution model with time-dependent diffusion-migration transport was developed for simulating localized corrosion. The mathematical formulation of this physical phenomenon and its computational formulation using a numerical scheme were described in detail to guide corrosion specialists. From a mathematical standpoint, the problem can be identified as a Stefan problem involving a convection-diffusion system. The model was applied to iron with the Butler-Volmer formula as the dissolution law, and both crevice and pit configurations were simulated. The results for crevice corrosion are qualitatively in good agreement with the predictions from the stationary crevice model made by Sharland, but our time-dependent model provides more
realistic results for concentration profiles and crevice depth. It was found that the dissolution rate increases with the chloride concentration and the metal potential, and that migration plays a significant role in the transport of species both in the crevice or the pit configuration. The evolution of the repassivation potential with respect to the pit depth was calculated and shown to decrease drastically and then more slowly, in good agreement with experimental results.

## Appendix A. ALE formulation

In this appendix, we present an Arbitrary Lagrangian Eulerian (ALE) formulation of the system (3.1),(2.22),(2.23) in order to handle the displacement of the free boundary (2.23). Let $t' \ge 0$ be fixed. The ALE formulation consists in transforming the current domain $[0, x_d(t)]$ for $t \ge 0$, onto the fixed interval $[0, x_d(t')]$ and to write down all the equations in this fixed interval. To do so, we define the ALE mapping (see [12, Chapter 5]), for $t \ge 0$:

$$ \psi(\cdot, t) : [0, x_d(t)] \rightarrow [0, x_d(t')] \\ x \mapsto \psi(x, t) = \frac{x_d(t')}{x_d(t)}x. \qquad (A.1) $$

The function $\psi$ is smooth with respect to $t$ and linear in $x$. This particular choice of $\psi$ is usually referred to as the Variable Space Grid (VSG) method [20, 16]. This mapping preserves the uniform discretization of the moving space interval $[0, x_d(t)]$ during the time evolution. In the system (3.1),(2.22),(2.23), we perform the change of unknowns

$$ \tilde{\phi}(\tilde{x}, t) = \phi(x, t) \quad \text{and} \quad \tilde{C}^j(\tilde{x}, t) = C_j(x, t) \quad \text{with} \quad \tilde{x} = \psi(x, t) \quad \text{for } t \ge 0, x \in [0, x_d(t)]. \qquad (A.2) $$

The new functions $\tilde{\phi}(\tilde{x}, t)$ and $\tilde{C}^j(\tilde{x}, t)$ together with the (unchanged) function $x_d$ satisfy the system:

**System satisfied by the concentrations $\{\tilde{C}^j\}_{j=1,...,M}$ and $\tilde{\phi}$**

$$ \frac{\partial \tilde{C}^j}{\partial t} = \mathcal{F}_1 \frac{\partial}{\partial \tilde{x}} \left( \hat{D}_j \left( \frac{\partial \tilde{C}^j}{\partial \tilde{x}} + z_j \tilde{C}^j \frac{\partial \tilde{\phi}}{\partial \tilde{x}} \right) \right) + \mathcal{F}_2 \frac{\partial \tilde{C}^j}{\partial \tilde{x}}, \quad \text{for } t > 0, 0 < \tilde{x} < x_d(t'), j \in [1, M], \qquad (A.3a) $$

$$ \sqrt{\mathcal{F}_1} \cdot \left( \hat{D}_1 \frac{\partial \tilde{C}^1}{\partial \tilde{x}} + z_1 \hat{D}_1 \tilde{C}^1 \frac{\partial \tilde{\phi}}{\partial \tilde{x}} \right) = (1 - \hat{\Omega}_{Fe} \tilde{C}^1) \hat{f}_{Fe^{2+}}(\tilde{\phi}, \tilde{C}^1), \quad \text{at } \tilde{x} = x_d(t') \quad (j = 1), \qquad (A.3b) $$

$$ \sqrt{\mathcal{F}_1} \cdot \left( \hat{D}_j \frac{\partial \tilde{C}^j}{\partial \tilde{x}} + z_j \hat{D}_j \tilde{C}^j \frac{\partial \tilde{\phi}}{\partial \tilde{x}} \right) = -\hat{\Omega}_{Fe} \tilde{C}^j \hat{f}_{Fe^{2+}}(\tilde{\phi}, \tilde{C}^1), \quad \text{at } \tilde{x} = x_d(t'), j \in [2, M], \qquad (A.3c) $$

$$ \tilde{C}^j(0, t) = \hat{\beta}_j, \quad t > 0, j \in [1, M], \qquad (A.3d) $$

$$ \tilde{C}^j(\tilde{x}, 0) = \tilde{C}^{j,0}(\tilde{x}), \quad \tilde{x} \in (0, x_d(t')), j \in [1, M], \qquad (A.3e) $$

$$ \sum_{j=1}^{M} z_j \tilde{C}^j(\tilde{x}, t) = 0, \quad t > 0, \quad 0 \le \tilde{x} \le x_d(t'), \qquad (A.4a) $$

$$ \tilde{\phi}(0, t) = 0, \quad t > 0, \qquad (A.4b) $$

along with the free boundary condition

$$ x_d'(t) = \hat{\Omega}_{Fe} \hat{f}_{Fe^{2+}}(\tilde{\phi}, \tilde{C}^1) \quad \text{at } \tilde{x} = x_d(t'), \qquad (A.5) $$
realistic results for concentration profiles and crevice depth. It was found that the dissolution rate increases with the chloride concentration and the metal potential, and that migration plays a significant role in the transport of species both in the crevice or the pit configuration. The evolution of the repassivation potential with respect to the pit depth was calculated and shown to decrease drastically and then more slowly, in good agreement with experimental results.

## Appendix A. ALE formulation

In this appendix, we present an Arbitrary Lagrangian Eulerian (ALE) formulation of the system (3.1),(2.22),(2.23) in order to handle the displacement of the free boundary (2.23). Let $t' \ge 0$ be fixed. The ALE formulation consists in transforming the current domain $[0, x_d(t)]$ for $t \ge 0$, onto the fixed interval $[0, x_d(t')]$ and to write down all the equations in this fixed interval. To do so, we define the ALE mapping (see [12, Chapter 5]), for $t \ge 0$:

$$ \psi(\cdot, t) : [0, x_d(t)] \rightarrow [0, x_d(t')] \\ x \mapsto \psi(x, t) = \frac{x_d(t')}{x_d(t)}x. \qquad (A.1) $$

The function $\psi$ is smooth with respect to $t$ and linear in $x$. This particular choice of $\psi$ is usually referred to as the Variable Space Grid (VSG) method [20, 16]. This mapping preserves the uniform discretization of the moving space interval $[0, x_d(t)]$ during the time evolution. In the system (3.1),(2.22),(2.23), we perform the change of unknowns

$$ \tilde{\phi}(\tilde{x}, t) = \phi(x, t) \quad \text{and} \quad \tilde{C}^j(\tilde{x}, t) = C_j(x, t) \quad \text{with} \quad \tilde{x} = \psi(x, t) \quad \text{for } t \ge 0, x \in [0, x_d(t)]. \qquad (A.2) $$

The new functions $\tilde{\phi}(\tilde{x}, t)$ and $\tilde{C}^j(\tilde{x}, t)$ together with the (unchanged) function $x_d$ satisfy the system:

System satisfied by the concentrations $\{\tilde{C}^j\}_{j=1,...,M}$ and $\tilde{\phi}$.

$$ \frac{\partial \tilde{C}^j}{\partial t} = \mathcal{F}_1 \frac{\partial}{\partial \tilde{x}} \left( \hat{D}_j \left( \frac{\partial \tilde{C}^j}{\partial \tilde{x}} + z_j \tilde{C}^j \frac{\partial \tilde{\phi}}{\partial \tilde{x}} \right) \right) + \mathcal{F}_2 \frac{\partial \tilde{C}^j}{\partial \tilde{x}}, \quad \text{for } t > 0, 0 < \tilde{x} < x_d(t'), j \in [1, M], \qquad (A.3a) $$

$$ \sqrt{\mathcal{F}_1} \cdot \left( \hat{D}_1 \frac{\partial \tilde{C}^1}{\partial \tilde{x}} + z_1 \hat{D}_1 \tilde{C}^1 \frac{\partial \tilde{\phi}}{\partial \tilde{x}} \right) = (1 - \hat{\Omega}_{Fe} \tilde{C}^1) \hat{f}_{Fe^{2+}}(\tilde{\phi}, \tilde{C}^1), \quad \text{at } \tilde{x} = x_d(t') \quad (j=1), \qquad (A.3b) $$

$$ \sqrt{\mathcal{F}_1} \cdot \left( \hat{D}_j \frac{\partial \tilde{C}^j}{\partial \tilde{x}} + z_j \hat{D}_j \tilde{C}^j \frac{\partial \tilde{\phi}}{\partial \tilde{x}} \right) = - \hat{\Omega}_{Fe} \tilde{C}^j \hat{f}_{Fe^{2+}}(\tilde{\phi}, \tilde{C}^1), \quad \text{at } \tilde{x} = x_d(t'), j \in [2, M], \qquad (A.3c) $$

$$ \tilde{C}^j(0, t) = \hat{\beta}_j, \quad t > 0, j \in [1, M], \qquad (A.3d) $$

$$ \tilde{C}^j(\tilde{x}, 0) = \tilde{C}^{j,0}(\tilde{x}), \quad \tilde{x} \in (0, x_d(t')), \quad j \in [1, M], \qquad (A.3e) $$

$$ \sum_{j=1}^{M} z_j \tilde{C}^j(\tilde{x}, t) = 0, \quad t > 0, \quad 0 \le \tilde{x} \le x_d(t'), \qquad (A.4a) $$

$$ \tilde{\phi}(0, t) = 0, \quad t > 0, \qquad (A.4b) $$

along with the free boundary condition

$$ x_d'(t) = \hat{\Omega}_{Fe} \hat{f}_{Fe^{2+}}(\tilde{\phi}, \tilde{C}^1) \quad \text{at } \tilde{x} = x_d(t'), \qquad (A.5) $$
where

$$ \mathcal{F}_1 = \mathcal{F}_1(t) = \left( \frac{x_d(t')}{x_d(t)} \right)^2, \quad \mathcal{F}_2 = \mathcal{F}_2(\tilde{x}, t) = \frac{x_d'(t)}{x_d(t)} \tilde{x}. \qquad (A.6) $$

and $\hat{f}_{\text{Fe}^{2+}}(\tilde{\phi}, \tilde{C}^1)$ is given in (2.17), namely

$$ \hat{f}_{\text{Fe}^{2+}}(\tilde{\phi}, \tilde{C}^1) = \hat{k}_a \exp(\hat{\phi}_m - \tilde{\phi}) - \hat{k}_c \frac{C_0}{C_{\text{ref}}} \tilde{C}^1 \exp(-(\hat{\phi}_m - \tilde{\phi})), \qquad (A.7) $$

where the constants $\hat{\phi}_m, \hat{k}_a$ and $\hat{k}_c$ in (A.7) are given in (2.18). The initial concentrations $\tilde{C}^{j,0}(\tilde{x}) = C_j^0(x)$ in (A.3e) are given functions. A numerical scheme for solving the ALE formulation (A.3)-(A.5) is described in the next Appendix B.

## Appendix B. Semi-implicit scheme for the ALE formulation

In this appendix, we describe the time and space approximations for the ALE formulation (A.3)-(A.5), leading to a time semi-implicit finite differences scheme. Let $\Delta t > 0$ be the discretization time step such that $t^n = n \Delta t$ be the discrete times for $n \in \mathbb{N}^*$. The time increment $\Delta t = t^n - t^{n-1}$ varies with $n$. We start with an initial pit bottom position

$$ x_d^0 = x_d(0) \text{ given,} \qquad (B.1) $$

and we introduce $(N + 1)$ equidistant grid points for discretizing the initial interval $[0, x_d^0]$ with:

$$ x_i^0 = i \Delta x \quad \text{for } i = 0, 1, \dots, N, \quad \Delta x = \frac{x_d^0}{N} \quad \text{and} \quad x_N^0 = x_d^0. \qquad (B.2) $$

The initial data are also given for the concentrations and the potential:

$$ \begin{gathered} C_i^{j,0} = \tilde{C}^j(x_i^0, 0) \text{ given, for } i = 0, \dots, N; \quad j = 1, \dots, M \qquad (B.3) \\ \phi_0^0 = 0 \quad \text{and} \quad \phi_N^0 = \tilde{\phi}(x_d^0, 0) \text{ given.} \end{gathered} $$

Let $(\{\tilde{C}^j\}_{j=1,...,M}, \tilde{\phi}, x_d)$ be the solution of (A.3)-(A.5) at time $t = t^n$ and suppose we know an approximation $x_d^n$ of the pit bottom position at time $t^n$ together with a discretization of the interval $[0, x_d^n]$ with $(N + 1)$ grid points $\tilde{x}_i$ ($i = 0, \dots, N$):

$$ x_d^n \approx x_d(t^n), \quad 0 = \tilde{x}_0 := x_0^n < \tilde{x}_1 := x_1^n < \dots < \tilde{x}_N := x_N^n = x_d^n. $$

In the ALE scheme, the number $(N + 1)$ of grid points remains constant at each time step $t = t^n$, $n \ge 0$, but the grid points $x_i^n$ ($i = 0, \dots, N$) of the interval $[0, x_d^n]$ are moving with the time iteration $n$. The SVG method we use (see [20, 16]) preserves the equidistance between grid points (see (B.8),(B.9)).

We also consider the approximations for the concentrations and the potential:

$$ C_i^{j,n} \approx \tilde{C}^j(\tilde{x}_i, t^n), \quad C_i^{j,n+1} \approx \tilde{C}^j(\tilde{x}_i, t^{n+1}) \quad \text{and} \quad \phi_i^n \approx \tilde{\phi}(\tilde{x}_i, t^n). \qquad (B.4) $$

for $j = 1, \dots, M$; $i = 0, \dots, N$.

Starting with the initial data (B.1)-(B.3), we compute $(\{C_i^{j,n+1}\}_{j=1,...,M}, \phi_i^{n+1}, x_d^{n+1})$ from $(\{C_i^{j,n}\}_{j=1,...,M}, \phi_N^n, x_d^n)$ for $n \ge 0$ and $i = 0, \dots, N$, according to the decoupling scheme presented
where

$$ \mathcal{F}_1 = \mathcal{F}_1(t) = \left( \frac{x_d(t')}{x_d(t)} \right)^2, \quad \mathcal{F}_2 = \mathcal{F}_2(\tilde{x}, t) = \frac{x'_d(t)}{x_d(t)} \tilde{x}. \qquad (A.6) $$

and $\hat{f}_{\text{Fe}^{2+}}(\tilde{\phi}, \tilde{C}^1)$ is given in (2.17), namely

$$ \hat{f}_{\text{Fe}^{2+}}(\tilde{\phi}, \tilde{C}^1) = \hat{k}_a \exp(\hat{\phi}_m - \tilde{\phi}) - \hat{k}_c \frac{C_0}{C_{\text{ref}}} \tilde{C}^1 \exp(-(\hat{\phi}_m - \tilde{\phi})), \qquad (A.7) $$

where the constants $\hat{\phi}_m, \hat{k}_a$ and $\hat{k}_c$ in (A.7) are given in (2.18). The initial concentrations $\tilde{C}^{j,0}(\tilde{x}) = C_j^0(x)$ in (A.3e) are given functions. A numerical scheme for solving the ALE formulation (A.3)-(A.5) is described in the next Appendix B.

## Appendix B. Semi-implicit scheme for the ALE formulation

In this appendix, we describe the time and space approximations for the ALE formulation (A.3)-(A.5), leading to a time semi-implicit finite differences scheme. Let $\Delta t > 0$ be the discretization time step such that $t^n = n \Delta t$ be the discrete times for $n \in \mathbb{N}^*$. The time increment $\Delta t = t^n - t^{n-1}$ varies with $n$. We start with an initial pit bottom position

$$ x_d^0 = x_d(0) \text{ given,} \qquad (B.1) $$

and we introduce $(N + 1)$ equidistant grid points for discretizing the initial interval $[0, x_d^0]$ with:

$$ x_i^0 = i \Delta x \quad \text{for } i = 0, 1, \dots, N, \quad \Delta x = \frac{x_d^0}{N} \quad \text{and} \quad x_N^0 = x_d^0. \qquad (B.2) $$

The initial data are also given for the concentrations and the potential:

$$ \begin{gathered} C_i^{j,0} = \tilde{C}^j(x_i^0, 0) \text{ given, for } i = 0, \dots, N; \quad j = 1, \dots, M \qquad (B.3) \\ \phi_0^0 = 0 \quad \text{and} \quad \phi_N^0 = \tilde{\phi}(x_d^0, 0) \text{ given.} \end{gathered} $$

Let $(\{\tilde{C}^j\}_{j=1,...,M}, \tilde{\phi}, x_d)$ be the solution of (A.3)-(A.5) at time $t = t^n$ and suppose we know an approximation $x_d^n$ of the pit bottom position at time $t^n$ together with a discretization of the interval $[0, x_d^n]$ with $(N + 1)$ grid points $\tilde{x}_i$ ($i = 0, \dots, N$):

$$ x_d^n \approx x_d(t^n), \quad 0 = \tilde{x}_0 := x_0^n < \tilde{x}_1 := x_1^n < \dots < \tilde{x}_N := x_N^n = x_d^n. $$

In the ALE scheme, the number $(N + 1)$ of grid points remains constant at each time step $t = t^n$, $n \ge 0$, but the grid points $x_i^n$ ($i = 0, \dots, N$) of the interval $[0, x_d^n]$ are moving with the time iteration $n$. The SVG method we use (see [20, 16]) preserves the equidistance between grid points (see (B.8), (B.9)).

We also consider the approximations for the concentrations and the potential:

$$ C_i^{j,n} \approx \tilde{C}^j(\tilde{x}_i, t^n), \quad C_i^{j,n+1} \approx \tilde{C}^j(\tilde{x}_i, t^{n+1}) \quad \text{and} \quad \phi_i^n \approx \tilde{\phi}(\tilde{x}_i, t^n). \qquad (B.4) $$

for $j = 1, \dots, M$; $i = 0, \dots, N$.

Starting with the initial data (B.1)-(B.3), we compute $(\{C_i^{j,n+1}\}_{j=1,...,M}, \phi_i^{n+1}, x_d^{n+1})$ from $(\{C_i^{j,n}\}_{j=1,...,M}, \phi_N^n, x_d^n)$ for $n \ge 0$ and $i = 0, \dots, N$, according to the decoupling scheme presented
below for the ALE system (A.3)–(A.5). The procedure to determine the initial potential $\phi_i^0$ for all $i = 1, \dots, N-1$ from the initial data (B.1)–(B.3), is described in subsection 3.2. We define the forward and backward differencing operators by:

$$ (\mathcal{D}_{\Delta x}^{+} u)_i = \frac{u_{i+1} - u_i}{\Delta x}, \quad (\mathcal{D}_{\Delta x}^{-} u)_i = \frac{u_i - u_{i-1}}{\Delta x}. \qquad (B.5) $$

* **Step 1: Interface displacement** (A.5). We start by computing the new position $x_d^{n+1}$ of the interface using the Butler-Volmer formula

$$ x_d^{n+1} = x_d^n + v^n \Delta t \quad \text{with} \quad v^n = \hat{\Omega}_{Fe} \hat{f}_{Fe^{2+}}(\phi_N^n, C_N^{1,n}) \quad \text{at} \quad x = x_d^n = x_N^n, \qquad (B.6) $$

i.e. with $v^n \simeq x'_d(t^n)$ (forward Euler), where $\hat{f}_{Fe^{2+}}(\phi_N^n, C_N^{1,n})$ is given by (see (A.7))

$$ \hat{f}_{Fe^{2+}}(\phi_N^n, C_N^{1,n}) = \hat{k}_a \exp(\hat{\phi}_m - \phi_N^n) - \hat{k}_c \frac{C_0}{C_{\text{ref}}} C_N^{1,n} \exp(-(\hat{\phi}_m - \phi_N^n)). \qquad (B.7) $$

We choose $\Delta t$ small enough so that $v^n \Delta t$ does not exceed $\Delta x$.

* **Step 2: Computation of the new mesh of $[0, x_d^{n+1}]$.** Taking $t = t^n$ and $t' = t^{n+1}$ in (A.1) leads to introduce

$$ \begin{array}{rcl} \psi^{n+1} : [0, x_d^n] & \to & [0, x_d^{n+1}] \\ x & \mapsto & \tilde{x} = \frac{x_d^{n+1}}{x_d^n} x & \qquad (B.8) \end{array} $$

which maps $[0, x_d^n]$ onto $[0, x_d^{n+1}]$. The $(N+1)$ grid points of the interval $[0, x_d^{n+1}]$ are given by:

$$ x_i^{n+1} := \tilde{x}_i = \psi^{n+1}(x_i^n), \quad \text{for } i = 0, \dots, N. \qquad (B.9) $$

* **Step 3: Full implicit in time and space discretizations of system (A.3),(A.4): computation of $C_i^{j,n+1}$ and $\phi_i^{n+1}$ for $j = 1, \dots, M$; $i = 0, \dots, N$**.

We take $t' = t^{n+1}$ in (A.3) at time $t = t^{n+1}$. We solve the following problem for $C_i^{j,n+1}$ and $\phi_i^{n+1}$ for $i = 0, \dots, N$:

- The electrodiffusion equations (A.3a) are discretized as:

$$ \frac{C_i^{j,n+1} - C_i^{j,n}}{\Delta t} = - (\mathcal{D}_{\Delta x}^{-} \mathbf{J}^{j,n+1})_i + \frac{\hat{\Omega}_{Fe} \hat{f}_{Fe^{2+}}(\phi_N^{n+1}, C_N^{1,n+1})}{x_d^{n+1}} x_i^{n+1} (\mathcal{D}_{\Delta x}^{-} C^{j,n+1})_i, \qquad (B.10) $$

with

$$ \mathbf{J}_i^{j,n+1} = -\hat{D}_j \left( (\mathcal{D}_{\Delta x}^{+} C^{j,n+1})_i + z_j \left( \frac{C_{i+1}^{j,n+1} + C_i^{j,n+1}}{2} \right) (\mathcal{D}_{\Delta x}^{+} \phi^{n+1})_i \right) \qquad (B.11) $$

for $j = 1, \dots, M$ and $i = 1, \dots, N-1$

- Boundary condition at $x = x_N^{n+1} = x_d^{n+1}$ for the metallic cation (see (A.3b)):

$$ \mathbf{J}_N^{1,n+1} = -(1 - \hat{\Omega}_{Fe} C_N^{1,n+1}) \hat{f}_{Fe^{2+}}(\phi_N^{n+1}, C_N^{1,n+1}). \qquad (B.12) $$
below for the ALE system (A.3)–(A.5). The procedure to determine the initial potential $\phi_i^0$ for all $i = 1, \dots, N-1$ from the initial data (B.1)–(B.3), is described in subsection 3.2. We define the forward and backward differencing operators by:

$$ (\mathcal{D}_{\Delta x}^{+} u)_i = \frac{u_{i+1} - u_i}{\Delta x}, \quad (\mathcal{D}_{\Delta x}^{-} u)_i = \frac{u_i - u_{i-1}}{\Delta x}. \qquad (B.5) $$

* **Step 1: Interface displacement** (A.5). We start by computing the new position $x_d^{n+1}$ of the interface using the Butler-Volmer formula

$$ x_d^{n+1} = x_d^n + v^n \Delta t \quad \text{with} \quad v^n = \hat{\Omega}_{Fe} \hat{f}_{Fe^{2+}}(\phi_N^n, C_N^{1,n}) \quad \text{at} \quad x = x_d^n = x_N^n, \qquad (B.6) $$

i.e. with $v^n \simeq x'_d(t^n)$ (forward Euler), where $\hat{f}_{Fe^{2+}}(\phi_N^n, C_N^{1,n})$ is given by (see (A.7))

$$ \hat{f}_{Fe^{2+}}(\phi_N^n, C_N^{1,n}) = \hat{k}_a \exp(\hat{\phi}_m - \phi_N^n) - \hat{k}_c \frac{C_0}{C_{\text{ref}}} C_N^{1,n} \exp(-(\hat{\phi}_m - \phi_N^n)). \qquad (B.7) $$

We choose $\Delta t$ small enough so that $v^n \Delta t$ does not exceed $\Delta x$.

* **Step 2: Computation of the new mesh of $[0, x_d^{n+1}]$.** Taking $t = t^n$ and $t' = t^{n+1}$ in (A.1) leads to introduce

$$ \begin{array}{rcl} \psi^{n+1} : [0, x_d^n] & \to & [0, x_d^{n+1}] \\ x & \mapsto & \tilde{x} = \frac{x_d^{n+1}}{x_d^n} x & \qquad (B.8) \end{array} $$

which maps $[0, x_d^n]$ onto $[0, x_d^{n+1}]$. The $(N+1)$ grid points of the interval $[0, x_d^{n+1}]$ are given by:

$$ x_i^{n+1} := \tilde{x}_i = \psi^{n+1}(x_i^n), \quad \text{for } i = 0, \dots, N. \qquad (B.9) $$

* **Step 3: Full implicit in time and space discretizations of system (A.3),(A.4): computation of $C_i^{j,n+1}$ and $\phi_i^{n+1}$ for $j = 1, \dots, M$; $i = 0, \dots, N$**.

We take $t' = t^{n+1}$ in (A.3) at time $t = t^{n+1}$. We solve the following problem for $C_i^{j,n+1}$ and $\phi_i^{n+1}$ for $i = 0, \dots, N$:

- The electrodiffusion equations (A.3a) are discretized as:

$$ \frac{C_i^{j,n+1} - C_i^{j,n}}{\Delta t} = - (\mathcal{D}_{\Delta x}^{-} \mathbf{J}^{j,n+1})_i + \frac{\hat{\Omega}_{Fe} \hat{f}_{Fe^{2+}}(\phi_N^{n+1}, C_N^{1,n+1})}{x_d^{n+1}} x_i^{n+1} (\mathcal{D}_{\Delta x}^{-} C^{j,n+1})_i, \qquad (B.10) $$

with

$$ \mathbf{J}_i^{j,n+1} = -\hat{D}_j \left( (\mathcal{D}_{\Delta x}^{+} C^{j,n+1})_i + z_j \left( \frac{C_{i+1}^{j,n+1} + C_i^{j,n+1}}{2} \right) (\mathcal{D}_{\Delta x}^{+} \phi^{n+1})_i \right) \qquad (B.11) $$

for $j = 1, \dots, M$ and $i = 1, \dots, N-1$

- Boundary condition at $x = x_N^{n+1} = x_d^{n+1}$ for the metallic cation (see (A.3b)):

$$ \mathbf{J}_N^{1,n+1} = -(1 - \hat{\Omega}_{Fe} C_N^{1,n+1}) \hat{f}_{Fe^{2+}}(\phi_N^{n+1}, C_N^{1,n+1}). \qquad (B.12) $$
- Boundary conditions at $x = x_N^{n+1} = x_d^{n+1}$ for the other ions (see (A.3c)):

$$ \mathbf{J}_N^{j,n+1} = \hat{\Omega}_{Fe} C_N^{j,n+1} \hat{f}_{Fe^{2+}}(\phi_N^{n+1}, C_N^{1,n+1}), \quad \text{for } j = 2, \dots, M \quad (B.13) $$

where

$$ \mathbf{J}_N^{j,n+1} = -\hat{D}_j \left( \frac{3 C_N^{j,n+1} - 4 C_{N-1}^{j,n+1} + C_{N-2}^{j,n+1}}{2\Delta x} + z_j C_N^{j,n+1} \left( \frac{3\phi_N^{n+1} - 4\phi_{N-1}^{n+1} + \phi_{N-2}^{n+1}}{2\Delta x} \right) \right). \quad (B.14) $$

for $j = 1, \dots, M$.

- Boundary conditions at $x = x_0^{n+1} = 0$ (see (A.3d), (A.4b)):

$$ C_0^{j,n+1} = \hat{\beta}_j \text{ for } j = 1, \dots, M \text{ such that } \sum_{j=1}^{M} z_j C_0^{j,n+1} = 0, \phi_0^{n+1} = 0. \quad (B.15) $$

- Electroneutrality condition (see (A.4a)):

$$ \sum_{j=1}^{M} z_j C_i^{j,n+1} = 0 \text{ for } i = 1, \dots, N. \quad (B.16) $$

**Remark 1.** The convective term in (B.10) is treated with an upwind scheme involving the discrete operator $D^-_{\Delta x}$ since the velocity flux $\hat{f}_{Fe^{2+}}$ is expected to be positive.

**Remark 2.** Due to (B.8), we have that

$$ \begin{aligned} C_i^{j,n} &\approx \tilde{C}^j(\tilde{x}_i, t^n) = \tilde{C}^j(\psi^{n+1}(x_i^n), t^n) = C_j(x_i^n, t^n), \\ \phi_i^n &\approx \tilde{\phi}(\tilde{x}_i, t^n) = \tilde{\phi}(\psi^{n+1}(x_i^n), t^n) = \phi(x_i^n, t^n), \end{aligned} \quad (B.17) $$

Hence we see that $C_i^{j,n}$ and $\phi_i^n$ are also approximations of $C_j(x_i^n, t^n)$ and $\phi(x_i^n, t^n)$ respectively, where $C_j$ and $\phi$ are the solutions of the original system (3.1), (2.22), (2.23).

With $m = (M + 1) \times (N + 1)$, the nonlinear problem (B.10)–(B.16) forms a nonlinear system of $m$ equations with $m$ unknowns $S = (\{C_i^{j,n+1}\}_{j=1,...,M}, \phi_i^{n+1}$ for $i = 0, \dots, N)$. This problem can be written as $\mathcal{F}(S) = 0$ and we use a Newton method to numerically solve it, with $S^0 = (\{C_i^{j,n}\}_{j=1,...,M}, \phi_i^n$; for $i = 0, \dots, N)$ as the initial guess. We refer to [12] for more details.

## Appendix C. The diffusion-migration-reaction model for five ions

We present an extended model with additional reaction terms by considering the hydrolysis reactions of iron cations. This reaction induces the formation of $H^+$ cations which changes the acidity of the solution, a criterium often proposed in literature (see e.g. [27]) for the pit stability. The dimensional Nernst-Planck equation introduced in (2.1a) becomes:

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right) + \mathcal{R}_i(C_1, \dots, C_M), \quad j = 1, \dots, M. \quad (C.1) $$
- Boundary conditions at $x = x_N^{n+1} = x_d^{n+1}$ for the other ions (see (A.3c)):

$$ \mathbf{J}_N^{j,n+1} = \hat{\Omega}_{Fe} C_N^{j,n+1} \hat{f}_{Fe^{2+}} (\phi_N^{n+1}, C_N^{1,n+1}), \quad \text{for } j = 2, \dots, M \quad (B.13) $$

where

$$ \mathbf{J}_N^{j,n+1} = -\hat{D}_j \left( \frac{3 C_N^{j,n+1} - 4 C_{N-1}^{j,n+1} + C_{N-2}^{j,n+1}}{2\Delta x} + z_j C_N^{j,n+1} \left( \frac{3\phi_N^{n+1} - 4\phi_{N-1}^{n+1} + \phi_{N-2}^{n+1}}{2\Delta x} \right) \right). \quad (B.14) $$

for $j = 1, \dots, M$.

- Boundary conditions at $x = x_0^{n+1} = 0$ (see (A.3d), (A.4b)):

$$ C_0^{j,n+1} = \hat{\beta}_j \text{ for } j = 1, \dots, M \text{ such that } \sum_{j=1}^{M} z_j C_0^{j,n+1} = 0, \phi_0^{n+1} = 0. \quad (B.15) $$

- Electroneutrality condition (see (A.4a)):

$$ \sum_{j=1}^{M} z_j C_i^{j,n+1} = 0 \text{ for } i = 1, \dots, N. \quad (B.16) $$

**Remark 1.** The convective term in (B.10) is treated with an upwind scheme involving the discrete operator $\mathcal{D}_{\Delta x}^{-}$ since the velocity flux $\hat{f}_{Fe^{2+}}$ is expected to be positive.

**Remark 2.** Due to (B.8), we have that

$$ \begin{aligned} C_i^{j,n} &\approx \tilde{C}^j(\tilde{x}_i, t^n) = \tilde{C}^j(\psi^{n+1}(x_i^n), t^n) = C_j(x_i^n, t^n), \\ \phi_i^n &\approx \tilde{\phi}(\tilde{x}_i, t^n) = \tilde{\phi}(\psi^{n+1}(x_i^n), t^n) = \phi(x_i^n, t^n), \end{aligned} \quad (B.17) $$

Hence we see that $C_i^{j,n}$ and $\phi_i^n$ are also approximations of $C_j(x_i^n, t^n)$ and $\phi(x_i^n, t^n)$ respectively, where $C_j$ and $\phi$ are the solutions of the original system (3.1), (2.22), (2.23).

With $m = (M + 1) \times (N + 1)$, the nonlinear problem (B.10)–(B.16) forms a nonlinear system of $m$ equations with $m$ unknowns $S = (\{C_i^{j,n+1}\}_{j=1,...,M}, \phi_i^{n+1}$ for $i = 0, \dots, N)$. This problem can be written as $\mathcal{F}(S) = 0$ and we use a Newton method to numerically solve it, with $S^0 = (\{C_i^{j,n}\}_{j=1,...,M}, \phi_i^n$; for $i = 0, \dots, N)$ as the initial guess. We refer to [12] for more details.

## Appendix C. The diffusion-migration-reaction model for five ions

We present an extended model with additional reaction terms by considering the hydrolysis reactions of iron cations. This reaction induces the formation of $H^+$ cations which changes the acidity of the solution, a criterium often proposed in literature (see e.g. [27]) for the pit stability. The dimensional Nernst-Planck equation introduced in (2.1a) becomes:

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right) + \mathcal{R}_i(C_1, \dots, C_M), \quad j = 1, \dots, M. \quad (C.1) $$
where $R_i(C_1, \dots, C_M)$ stands for the reaction terms for the concentrations involved in the chemical reaction under consideration. In this appendix, we describe a model with five ions including diffusion, migration and reaction terms by considering the following chemical reaction:

$$ \mathrm{Fe}^{2+} + \mathrm{H}_2\mathrm{O} \underset{k_b}{\overset{k_f}{\rightleftharpoons}} \mathrm{Fe}(\mathrm{OH})^{+} + \mathrm{H}^{+}. \qquad (\mathrm{C.2}) $$

where $k_f$ and $k_b$ are the forward and backward reaction rate constants. For simplicity, we set

$$ C_1 := C_{Fe^{2+}}, \quad C_2 := C_{Na^{+}}, \quad C_3 := C_{Cl^{-}}, \quad C_4 := C_{H^{+}}, \quad C_5 := C_{Fe(OH)^{+}}. \qquad (\mathrm{C.3}) $$

When equilibrium is reached, the forward and backward rates in (C.2) become equal, and the amount of reactants and products obeys the fundamental law of mass action [28]:

$$ K_1 = \frac{k_f}{k_b} = \frac{\frac{C_4}{C_{ref}} \frac{C_5}{C_{ref}}}{\frac{C_1}{C_{ref}} \frac{C_{H_2O}}{C_{ref}}} \approx 10^{-9.316} \qquad (\mathrm{C.4}) $$

where $C_{H_2O} = C_{ref} = 1 \mathrm{mol/L}$ and $K_1$ is called the equilibrium constant of the chemical reaction (C.2). In order to take in account this chemical reaction, we have to modify the system (2.12)–(2.14),(2.7) by replacing the equations (2.13a) for the concentrations (without reaction terms) by the following set of equations including reaction terms:

$$ \frac{\partial C_1}{\partial t} = \frac{\partial}{\partial x} \left( D_1 \left( \frac{\partial C_1}{\partial x} + \frac{z_1 F}{RT} C_1 \frac{\partial \phi}{\partial x} \right) \right) - \frac{k_f}{C_{ref}} \left( C_1 - \frac{1}{K_1} \frac{C_4 C_5}{C_{ref}} \right), \qquad (\mathrm{C.5a}) $$

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right) \quad \text{for } j = 2, 3 \qquad (\mathrm{C.5b}) $$

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right) + \frac{k_f}{C_{ref}} \left( C_1 - \frac{1}{K_1} \frac{C_4 C_5}{C_{ref}} \right) \quad \text{for } j = 4, 5 \qquad (\mathrm{C.5c}) $$

The equilibrium state of the reaction (C.2) is achieved when $k_f \to +\infty$, leading to the equilibrium relation

$$ C_1 - \frac{1}{K_1} \frac{C_4 C_5}{C_{ref}} = 0. \qquad (\mathrm{C.6}) $$

From the numerical view point, we choose a very large value for the reaction rate constant $k_f$ to ensure that the equilibrium relation (C.6) is approximatively satisfies everywhere and at any time and also to ensure that the Newton method converges. For a value of $k_f \gg 10^5 \mathrm{mol.m}^{-3}.\mathrm{s}^{-1}$, we need to choose a very small time step, and as a result computation times are then very large. In the numerical simulations, we choose $k_f = 10^5 \mathrm{mol.m}^{-3}.\mathrm{s}^{-1}$ which seems to be a good compromise between satisfying the equilibrium relation (C.6) and the computational time cost. For more details and further numerical investigations we refer to [12, §6.6].

## Acknowledgements

The authors would like to thank Professor Yoichiro Mori from the University of Pennsylvania for stimulating discussions and useful comments. Moreover, they are very grateful to the referee of this paper who has helped us to improve the manuscript. This project was funded by CEA (COSTO project). Prof. H. Matano was also financially supported by KAKENHI 21H00995.
where $R_i(C_1, \dots, C_M)$ stands for the reaction terms for the concentrations involved in the chemical reaction under consideration. In this appendix, we describe a model with five ions including diffusion, migration and reaction terms by considering the following chemical reaction:

$$ \mathrm{Fe}^{2+} + \mathrm{H}_2\mathrm{O} \underset{k_b}{\overset{k_f}{\rightleftharpoons}} \mathrm{Fe}(\mathrm{OH})^{+} + \mathrm{H}^{+}. \qquad (\mathrm{C.2}) $$

where $k_f$ and $k_b$ are the forward and backward reaction rate constants. For simplicity, we set

$$ C_1 := C_{Fe^{2+}}, \quad C_2 := C_{Na^{+}}, \quad C_3 := C_{Cl^{-}}, \quad C_4 := C_{H^{+}}, \quad C_5 := C_{Fe(OH)^{+}}. \qquad (\mathrm{C.3}) $$

When equilibrium is reached, the forward and backward rates in (C.2) become equal, and the amount of reactants and products obeys the fundamental law of mass action [28]:

$$ K_1 = \frac{k_f}{k_b} = \frac{\frac{C_4}{C_{ref}} \frac{C_5}{C_{ref}}}{\frac{C_1}{C_{ref}} \frac{C_{H_2O}}{C_{ref}}} \approx 10^{-9.316} \qquad (\mathrm{C.4}) $$

where $C_{H_2O} = C_{ref} = 1 \mathrm{mol/L}$ and $K_1$ is called the equilibrium constant of the chemical reaction (C.2). In order to take in account this chemical reaction, we have to modify the system (2.12)–(2.14),(2.7) by replacing the equations (2.13a) for the concentrations (without reaction terms) by the following set of equations including reaction terms:

$$ \frac{\partial C_1}{\partial t} = \frac{\partial}{\partial x} \left( D_1 \left( \frac{\partial C_1}{\partial x} + \frac{z_1 F}{RT} C_1 \frac{\partial \phi}{\partial x} \right) \right) - \frac{k_f}{C_{ref}} \left( C_1 - \frac{1}{K_1} \frac{C_4 C_5}{C_{ref}} \right), \qquad (\mathrm{C.5a}) $$

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right) \quad \text{for } j = 2, 3 \qquad (\mathrm{C.5b}) $$

$$ \frac{\partial C_j}{\partial t} = \frac{\partial}{\partial x} \left( D_j \left( \frac{\partial C_j}{\partial x} + \frac{z_j F}{RT} C_j \frac{\partial \phi}{\partial x} \right) \right) + \frac{k_f}{C_{ref}} \left( C_1 - \frac{1}{K_1} \frac{C_4 C_5}{C_{ref}} \right) \quad \text{for } j = 4, 5 \qquad (\mathrm{C.5c}) $$

The equilibrium state of the reaction (C.2) is achieved when $k_f \to +\infty$, leading to the equilibrium relation

$$ C_1 - \frac{1}{K_1} \frac{C_4 C_5}{C_{ref}} = 0. \qquad (\mathrm{C.6}) $$

From the numerical view point, we choose a very large value for the reaction rate constant $k_f$ to ensure that the equilibrium relation (C.6) is approximatively satisfies everywhere and at any time and also to ensure that the Newton method converges. For a value of $k_f \gg 10^5 \mathrm{mol.m}^{-3}.\mathrm{s}^{-1}$, we need to choose a very small time step, and as a result computation times are then very large. In the numerical simulations, we choose $k_f = 10^5 \mathrm{mol.m}^{-3}.\mathrm{s}^{-1}$ which seems to be a good compromise between satisfying the equilibrium relation (C.6) and the computational time cost. For more details and further numerical investigations we refer to [12, §6.6].

## Acknowledgements

The authors would like to thank Professor Yoichiro Mori from the University of Pennsylvania for stimulating discussions and useful comments. Moreover, they are very grateful to the referee of this paper who has helped us to improve the manuscript. This project was funded by CEA (COSTO project). Prof. H. Matano was also financially supported by KAKENHI 21H00995.
References

[1] A. Turnbull. Review of modelling of pit propagation kinetics. *British Corrosion Journal*, 28(4):297-308, 1993.

[2] Sun, W., Wang, L., Wu, T. and Liu, G. An arbitrary Lagrangian-Eulerian model for modelling the time-dependent evolution of crevice corrosion. *Corrosion Science*, 78, 233-243, 2014.

[3] Xiao, J. and Chaudhuri, S. Predictive modeling of localized corrosion: An application to aluminum alloys. *Electrochimica Acta*, 56(16), 5630-5641, 2011.

[4] Ansari, T. Q., Xiao, Z., Hu, S., Li, Y., Luo, J. L. and Shi, S. Q. Phase-field model of pitting corrosion kinetics in metallic materials. *Npj Computational Materials*, 2018 4:1, 4(1), 1-9, 2018.

[5] Sun, X., Srinivasan, J., Kelly, R. G. and Duddu, R. Numerical investigation of critical electro-chemical factors for pitting corrosion using a multi-species reactive transport model. *Corrosion Science*, 179, 109130, 2021.

[6] S. M. Sharland and P. W. Tasker. A mathematical model of crevice and pitting corrosion—i. the physical model. *Corrosion Science*, 28(6):603–620, 1988.

[7] S. M. Sharland. A mathematical model of crevice and pitting corrosion—ii. the mathematical solution. *Corrosion Science*, 28(6):621–630, 1988.

[8] Brackman, M. D., Clemons, C. B., Golovaty, D., Kreider, K. L., Wilder, J., Young, G. W., Payer, J. and Lillard, R. S. Modeling and Simulation of Damage Evolution during Crevice Corrosion. *Journal of The Electrochemical Society*, 161(5), C237-C245, 2014.

[9] S. Scheiner and C. Hellmich. Stable pitting corrosion of stainless steel as diffusion-controlled dissolution process with a sharp moving electrode boundary. *Corrosion science*, 49(2):319-346, 2007.

[10] S. Scheiner and C. Hellmich. Finite volume model for diffusion-and activation-controlled pitting corrosion of stainless steel. *Computer Methods in Applied Mechanics and Engineering*, 198(37-40):2898-2910, 2009.

[11] C. V. Moraes and R. G. Kelly. A comparison of FEM results from the use of different governing equations in a galvanic cell part I: In the presence of a supporting electrolyte, *Electrochimica Acta*, Vol. 469, 143146, 2023.

[12] M. Bouguezzi. *Modeling and computer simulation of the propagation rate of pit corrosion*. PhD thesis, Paris-Saclay University, 2021. NNT : 2021UPASM027. tel-03456588f; https://theses.hal.science/tel-03456588v1

[13] J. Newman and K. E. Thomas-Alyea. Electrochemical systems. *John Wiley & Sons*, 2012.

[14] J. Caldwell and Y-Y. Kwan. A brief review of several numerical methods for one-dimensional Stefan problems. *Thermal Science*, 13(2):61-72, 2009.

[15] E. Javierre, C. Vuik, F.J. Vermolen, and S. Van der Zwaag. A comparison of numerical models for one-dimensional stefan problems. *Journal of Computational and Applied Mathematics*, 192(2):445-459, 2006.
References

[1] A. Turnbull. Review of modelling of pit propagation kinetics. *British Corrosion Journal*, 28(4):297-308, 1993.

[2] Sun, W., Wang, L., Wu, T. and Liu, G. An arbitrary Lagrangian-Eulerian model for modelling the time-dependent evolution of crevice corrosion. *Corrosion Science*, 78, 233-243, 2014.

[3] Xiao, J. and Chaudhuri, S. Predictive modeling of localized corrosion: An application to aluminum alloys. *Electrochimica Acta*, 56(16), 5630-5641, 2011.

[4] Ansari, T. Q., Xiao, Z., Hu, S., Li, Y., Luo, J. L. and Shi, S. Q. Phase-field model of pitting corrosion kinetics in metallic materials. *Npj Computational Materials*, 2018 4:1, 4(1), 1-9, 2018.

[5] Sun, X., Srinivasan, J., Kelly, R. G. and Duddu, R. Numerical investigation of critical electro-chemical factors for pitting corrosion using a multi-species reactive transport model. *Corrosion Science*, 179, 109130, 2021.

[6] S. M. Sharland and P. W. Tasker. A mathematical model of crevice and pitting corrosion—i. the physical model. *Corrosion Science*, 28(6):603–620, 1988.

[7] S. M. Sharland. A mathematical model of crevice and pitting corrosion—ii. the mathematical solution. *Corrosion Science*, 28(6):621–630, 1988.

[8] Brackman, M. D., Clemons, C. B., Golovaty, D., Kreider, K. L., Wilder, J., Young, G. W., Payer, J. and Lillard, R. S. Modeling and Simulation of Damage Evolution during Crevice Corrosion. *Journal of The Electrochemical Society*, 161(5), C237-C245, 2014.

[9] S. Scheiner and C. Hellmich. Stable pitting corrosion of stainless steel as diffusion-controlled dissolution process with a sharp moving electrode boundary. *Corrosion science*, 49(2):319–346, 2007.

[10] S. Scheiner and C. Hellmich. Finite volume model for diffusion-and activation-controlled pitting corrosion of stainless steel. *Computer Methods in Applied Mechanics and Engineering*, 198(37-40):2898-2910, 2009.

[11] C. V. Moraes and R. G. Kelly. A comparison of FEM results from the use of different governing equations in a galvanic cell part I: In the presence of a supporting electrolyte, *Electrochimica Acta*, Vol. 469, 143146, 2023.

[12] M. Bouguezzi. *Modeling and computer simulation of the propagation rate of pit corrosion*. PhD thesis, Paris-Saclay University, 2021. NNT : 2021UPASM027. tel-03456588f; https://theses.hal.science/tel-03456588v1

[13] J. Newman and K. E. Thomas-Alyea. Electrochemical systems. *John Wiley & Sons*, 2012.

[14] J. Caldwell and Y-Y. Kwan. A brief review of several numerical methods for one-dimensional Stefan problems. *Thermal Science*, 13(2):61–72, 2009.

[15] E. Javierre, C. Vuik, F.J. Vermolen, and S. Van der Zwaag. A comparison of numerical models for one-dimensional stefan problems. *Journal of Computational and Applied Mathematics*, 192(2):445–459, 2006.
[16] S. Kutluay, A.R. Bahadir, A. Özdeş. The numerical solution of one-phase classical Stefan problem, J. Comput. Appl. Math., 81(1):135-144, 1997.

[17] G. Segal, K. Vuik, and F. Vermolen. A conserving discretization for the free boundary in a two-dimensional Stefan problem. *Journal of Computational Physics*, 141(1):1-21, 1998.

[18] J. Donea, A. Huerta, J-Ph. Ponthot, and A. Rodriguez-Ferran. Arbitrary lagrangian-eulerian methods, volume 1 of encyclopedia of computational mechanics, chapter 14. *John Wiley & Sons Ltd*, 3:1-25, 2004.

[19] T. J. R. Hughes, W. Liu and T. K. Zimmerman. Lagrangian-Eulerian finite element formulation for incompressible viscous flows, Comput. Methods Appl. Mech. Engrg. 29:329-349, 1981.

[20] W. D. Murray and F. Landis. Numerical and machine solutions of transient heat-conduction problems involving melting or freezing: part i—method of analysis and sample solutions. *Journal of Heat Transfer*, 81(2):106-112, 1959.

[21] S. M. Sharland. *A Theoretical Study of Crevice and Pitting Corrosion in Steels*. PhD thesis, University of London, 1988.

[22] J. L. Mousson, B. Vuillemin, R. Oltra, D. Crusset, G. Santarini, and P. Combrade. Modeling of the propagation of crevice corrosion. 2004.

[23] https://www.aqion.de/site/194

[24] N. J. Laycock and S. P. White. Computer simulation of single pit propagation in stainless steel under potentiostatic control. *Journal of the Electrochemical Society*, 148(7):B264-B275, 2001.

[25] J. Srinivasan, C. Liu, and R. G. Kelly. Geometric evolution of flux from a corroding one-dimensional pit and its implications on the evaluation of kinetic parameters for pit stability. *Journal of The Electrochemical Society*, 163(10):C694, 2016.

[26] H. S. Isaacs, J. H. Cho, M. L. Rivers, and S. R. Sutton. In situ x-ray microprobe study of salt layers during anodic dissolution of stainless steel in chloride solution. *Journal of the Electrochemical Society*, 142(4):1111, 1995.

[27] B. Malki, and B. Baroux. Modeling of metastable pitting: Towards a better understanding of the effect of alloying elements. *ECS Trans.*, 3(31) 273-284 (2007).

[28] HSC Thermochemistry R software version 10.3.4, Metso Outotec Finland O.
[16] S. Kutluay, A.R. Bahadir, A. Özdeş. The numerical solution of one-phase classical Stefan problem, J. Comput. Appl. Math., 81(1):135-144, 1997.

[17] G. Segal, K. Vuik, and F. Vermolen. A conserving discretization for the free boundary in a two-dimensional Stefan problem. *Journal of Computational Physics*, 141(1):1-21, 1998.

[18] J. Donea, A. Huerta, J-Ph. Ponthot, and A. Rodriguez-Ferran. Arbitrary lagrangian-eulerian methods, volume 1 of encyclopedia of computational mechanics, chapter 14. *John Wiley & Sons Ltd*, 3:1-25, 2004.

[19] T. J. R. Hughes, W. Liu and T. K. Zimmerman. Lagrangian-Eulerian finite element formulation for incompressible viscous flows, Comput. Methods Appl. Mech. Engrg. 29:329-349, 1981.

[20] W. D. Murray and F. Landis. Numerical and machine solutions of transient heat-conduction problems involving melting or freezing: part i—method of analysis and sample solutions. *Journal of Heat Transfer*, 81(2):106-112, 1959.

[21] S. M. Sharland. *A Theoretical Study of Crevice and Pitting Corrosion in Steels*. PhD thesis, University of London, 1988.

[22] J. L. Mousson, B. Vuillemin, R. Oltra, D. Crusset, G. Santarini, and P. Combrade. Modeling of the propagation of crevice corrosion. 2004.

[23] https://www.aqion.de/site/194

[24] N. J. Laycock and S. P. White. Computer simulation of single pit propagation in stainless steel under potentiostatic control. *Journal of the Electrochemical Society*, 148(7):B264-B275, 2001.

[25] J. Srinivasan, C. Liu, and R. G. Kelly. Geometric evolution of flux from a corroding one-dimensional pit and its implications on the evaluation of kinetic parameters for pit stability. *Journal of The Electrochemical Society*, 163(10):C694, 2016.

[26] H. S. Isaacs, J. H. Cho, M. L. Rivers, and S. R. Sutton. In situ x-ray microprobe study of salt layers during anodic dissolution of stainless steel in chloride solution. *Journal of the Electrochemical Society*, 142(4):1111, 1995.

[27] B. Malki, and B. Baroux. Modeling of metastable pitting: Towards a better understanding of the effect of alloying elements. *ECS Trans.*, 3(31) 273-284 (2007).

[28] HSC Thermochemistry R software version 10.3.4, Metso Outotec Finland O.