domain: phys math

text:
Meriem Bouguezzi

Jean-Francois Scheid

Danielle Hilhorst

Hiroshi Matano

Christian Bataillon

Florence Lequien

Fabien Rouillard

Jean-François Scheid

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

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

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 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.

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.

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 φ depend on a one space variable x, so that the system reduces to a one-dimensional problem. It is composed of the following domains :

Bulk solution

Metal

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, . . . , M ) at the position x and time t. The concentrations C j and the electrostatic potential φ = φ(x, t) satisfy the following Poisson-Nernst-Planck system [1] :

           ∂C j ∂t = ∂ ∂x D j ∂C j ∂x + z j F RT C j ∂φ ∂x , t > 0, 0 < x < x d (t), j = 1, . . . , M, (2.1a) - ∂ ∂x ∂φ ∂x = M j=1 z j F C j , t > 0, 0 < x < x d (t), (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 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 = -∂φ ∂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 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

M j=1 z j C j = 0. (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, . . . , M . With equation (2.2), this leads to

z 1 F ∂C j ∂t + • • • + z M F ∂C M ∂t = F ∂ ∂t (z 1 C 1 + • • • + z M C M ) = 0.
Consequently, we obtain from equation (2.1a),

M j=1 ∂ ∂x D j z j F ∂C j ∂x + z 2 j F 2 RT C j ∂φ ∂x = 0.
This equation can be rewritten as

∂ ∂x σ ∂φ ∂x + M j=1 ∂ ∂x z j F D j ∂C j ∂x = 0, t > 0, 0 < x < x d (t) where σ := M j=1 z 2 j F 2 RT D j C j , (2.3)
and σ is the electrical conductivity of the ionic solution.

Thus, in order to keep the electroneutrality of the solution, φ has to satisfy (2.3). In the special case where C j are spatially uniform (i.e. ∂C j ∂x = 0, j = 1, . . . , M ), equation (2.3) reduces to Laplace's equation ∂ 2 φ ∂x 2 = 0, hence φ is harmonic. However, φ 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).

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 φ, this is a free boundary problem. In other words, the unknowns are not just {C j } and φ, 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 :

x d (0) = x 0 d , C j (x, 0) = C 0 j (x) in [0, x 0 d (j = 1, . . . , M ), (2.4)
where x 0 d is the initial shape of the pit/crevice as shown in Figure 2.1, and C 0 j (x) (j = 1, . . . , M ) are functions satisfying

C 0 j ≥ 0 (j = 1, . . . , M ), z 1 C 0 1 + • • • + z M C 0 M = 0. (2.5)
We do not specify the initial value of φ since φ 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∈(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 2 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 ka ---- kc Fe 2+ + 2 e -.
To derive boundary conditions at the pit/crevice bottom, we need to distinguish Fe 2+ and other ions such as Na + and Cl -. In what follows, we let j = 1 represent Fe 2+ , that is,

C 1 := C Fe 2+ , z 1 := z Fe 2+ = +2.
Then the boundary conditions at x = x d (t) are given by

     J 1 = -f Fe 2+ + dx d (t) dt C 1 on x = x d (t), (2.6a)
J j = dx d (t) dt C j on x = x d (t) (j = 2, . . . , M ), (2.6b)
where 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 2+ . 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 2+ :

J j = -D j ∂C j ∂x + z j F RT C j ∂φ ∂x ,
f Fe 2+ (φ, C Fe 2+ ) = k a exp F RT (φ m -φ) -k c C Fe 2+ C ref exp - F RT (φ m -φ) , (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 2 s). C ref is a reference concentration which is usually set at C ref = 1 mol/L in ideal system, while φ m denotes the electric potential applied on the iron steel surface (which is assumed to be constant) and φ, 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 φ at the pit/crevice bottom x = x d (t). We can derive the boundary condition for φ 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, . . . , M ) and sum them up, and use (2.2). Then we obtain

σ ∂φ ∂x + M j=1 z j F D j ∂C j ∂x = z 1 F f Fe 2+ (φ, C 1 ) on x = x d (t), (2.8)
where σ

= M j=1 z 2 j F 2 RT D j C j .
Here C j , j = 1, . . . , M , are regarded as given functions. This is a nonlinear Robin boundary condition for φ.

Boundary conditions satisfied by {C j } j∈(1,••• ,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 = β j at x = 0, t > 0 (j = 1, 2, . . . , M ), (2.9)
where β j (j = 1, . . . , M ) are constants which satisfy the local electroneutrality condition

M j=1 z j β j = 0.
Boundary conditions satisfied by φ at the entrance x = 0 of the pit/crevice. Unlike the pit/crevice bottom, where we do not specify the boundary condition for φ, we need to specify the boundary condition for φ at x = 0 in order to determine φ 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]): φ = 0 at x = 0, t > 0.

(2.10)

Boundary condition at the moving interface x = x d (t).

The corrosion rate dx d (t) dt , namely the normal velocity of the boundary x d (t) in the outward direction, is proportional to the flux f Fe 2+ through the following relationship :

dx d (t) dt = x d (t) = Ω Fe f Fe 2+ (m/s), (2.11)
where Ω Fe is the molar volume of the solid iron, Ω Fe ≈ 7.1 cm 3 /mol = 7.1 • 10 -6 m 3 /mol.

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 φ and the concentrations {C j } j∈(1,••• ,M ) are given by: System satisfied by the potential φ

                   ∂ ∂x σ ∂φ ∂x + M j=1 ∂ ∂x z j F D j ∂C j ∂x = 0, t > 0, 0 < x < x d (t), (2.12a)
σ ∂φ ∂x + M j=1 z j F D j ∂C j ∂x = z 1 F f Fe 2+ (φ, C 1 ) at x = x d (t), (2.12b) φ(0, t) = 0, t > 0, (2.12c)
where σ is given in (2.3) and

f Fe 2+ is given in (2.7).
System satisfied by the concentrations

{C j } j=1,••• ,M                            ∂C j ∂t = ∂ ∂x D j ∂C j ∂x + z j F RT C j ∂φ ∂x , t > 0, 0 < x < x d (t), j = 1, . . . , M, (2.13a
)

D 1 ∂C 1 ∂x + z 1 F RT C 1 ∂φ ∂x = -x d (t)C 1 + f Fe 2+ (φ, C 1 ) at x = x d (t), (2.13b)
D j ∂C j ∂x + z j F RT C j ∂φ ∂x = -x d (t)C j at x = x d (t) (j = 2, . . . , M ), (2.13c) C j (0, t) = β j , t > 0 (j = 1, 2, . . . , M ), (2.13d)
C j (x, 0) = C 0 j (x), x ∈ [0, x 0 d ) (2.13e)
along with the free boundary condition

x d (t) = Ω Fe f Fe 2+ φ x d (t), t , C 1 x d (t), t at x = x d (t), x d (0) = x 0 d . (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 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 τ and L respectively, where τ = 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 x, t = L 2 D 0 t, C j = C 0 C j , D j = D 0 D j , φ = RT F φ.
(2.15)

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

System satisfied by the non-dimensional potential φ

                   ∂ ∂ x σ ∂ φ ∂ x + M j=1 ∂ ∂ x z j D j ∂ C j ∂ x = 0, t > 0, 0 < x < x d ( t), (2.16a) σ ∂ φ ∂ x + M j=1 z j D j ∂ C j ∂ x = z 1 f Fe 2+ ( φ, C 1 ) at x = x d ( t), C 1 = C Fe 2+ , (2.16b) φ(0, t) = 0, t > 0, (2.16c)
where σ := M j=1 z 2 j D j C j and the non-dimensional form of f Fe 2+ denoted by f Fe 2+ is given by

f Fe 2+ φ, C 1 = k a exp φ m -φ -k c C 1 exp -( φ m -φ) , (2.17)
where

φ m = F RT φ m , k a = L D 0 C 0 k a and k c = L D 0 C ref k c . (2.18)
System satisfied by the non-dimensional concentrations

{ C j } j=1,••• ,M                                      ∂ C j ∂ t = ∂ ∂ x D j ∂ C j ∂ x + z j C j ∂ φ ∂ x , t > 0, 0 < x < x d ( t), j = 1, . . . , M, (2.19a
)

D 1 ∂ C 1 ∂ x + z 1 C 1 ∂ φ ∂ x = -x d ( t) C 1 + f Fe 2+ φ, C 1 at x = x d ( t), (2.19b)
D j ∂ C j ∂ x + z j C j ∂ φ ∂ x = -x d ( t) C j at x = x d ( t) (j = 2, . . . , M ), (2.19c) C j (0, t) = β j , t > 0 (j = 1, 2, . . . , M ) where β j = β j C 0 , (2.19d)
C j (x, 0) = C 0 j (x), x ∈ [0, x 0 d ) (2.19e)
along with the non-dimensional free boundary condition

x d ( t) = Ω Fe f Fe 2+ φ x d ( t), t , C 1 x d ( t), t at x = x d ( t)
where

Ω Fe = C 0 Ω Fe , x d (0) = x 0 d . (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 • from variables C j , φ, x, t. Henceforth, the following non-dimensional system is considered from now on:

                   ∂ ∂x σ ∂φ ∂x + M j=1 ∂ ∂x z j D j ∂C j ∂x = 0, t > 0, 0 < x < x d (t), (2.21a)
σ ∂φ ∂x + M j=1 z j D j ∂C j ∂x = z 1 f Fe 2+ φ, C 1 at x = x d (t), (2.21b) φ(0, t) = 0, t > 0, (2.21c)                          ∂C j ∂t = ∂ ∂x D j ∂C j ∂x + z j C j ∂φ ∂x , t > 0, 0 < x < x d (t), (j = 1, • • • , M ) (2.22a) D 1 ∂C 1 ∂x + z 1 D 1 C 1 ∂φ ∂x = -x d (t)C 1 + f Fe 2+ (φ, C 1 ) at x = x d (t), (2.22b)
D j ∂C j ∂x + z j D j C j ∂φ ∂x = -x d (t)C j at x = x d (t) (j = 2, • • • , M ), (2.22c) C j (0, t) = β j , t > 0 (j = 1, • • • , M ), (2.22d) C j (x, 0) = C 0 j (x), x ∈ [0, x 0 d ), (2.22e)
along with the free boundary condition

x d (t) = Ω Fe f Fe 2+ φ x d (t), t , C 1 x d (t), t at x = x d (t), x d (0) = x 0 d , (2.23)
where

f Fe 2+ (φ, C 1 ) = k a exp φ m -φ -k c C 1 exp -( φ m -φ) .
(2.24)

3. Numerical approximation of the one-dimensional model

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 diffusionconvection 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:

     M j=1 z j C j (x, t) = 0, t > 0, 0 x x d (t), (3.1a) φ(0, t) = 0, t > 0. (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.

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 φ is described below. At t = 0, the only unknown in (2.21) is the potential

φ since C 1 , C 2 , • • • , C M , x d
are given functions at this stage. In fact, for all t 0, the potential φ can be determined as a function of C 1 , C 2 , • • • , C M and x d . Indeed, integrating (2.21a) by x shows that

σ ∂φ ∂x + M j=1 z j D j ∂C j ∂x (3.2)
is independent of x, hence, by (2.21b), the following equation is obtained

σ ∂φ ∂x + M j=1 z j D j ∂C j ∂x = z 1 f Fe 2+ P (t), C 1 (x d (t), t) for 0 ≤ x ≤ x d (t), t > 0 (3.3)
where

P (t) := φ x d (t), t . (3.4)
Then (3.3) can be rewritten as

∂φ ∂x = - M j=1 z j D j ∂C j ∂x M j=1 z 2 j D j C j + z 1 f Fe 2+ P (t), C 1 (x d (t), t) M j=1 z 2 j D j C j . (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 f Fe 2+ P (t), C 1 (x d (t), t) B(t), (3.6)
where

A(t) = x d (t) 0 M j=1 z j D j ∂C j ∂x M j=1 z 2 j D j C j dx, B(t) = x d (t) 0 dx M j=1 z 2 j D j C j = x d (t) 0 dx σ(x, t) . (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 0, the potential at the pit/crevice bottom, P (t) = φ(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∈{1,••• ,M } and x d ), the derivative ∂φ ∂x is determined on the entire interval 0 x x d (t) from (3.5). Thus, the whole initial potential 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.

φ 0 = φ( • , 0) in [0, x 0 d ] at t = 0 can be computed from the initial data x 0 d and C 0 j , (j = 1, • • • , M ).
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):

Fe 2+ + H 2 O kf ---- kb Fe (OH) + + H + . (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-migrationreaction 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 • 10 -9 m 2 .s -1 , L = 10 -3 m, C 0 = 10 -3 mol/L.
The diffusion coefficients of concerned species are given by (see [23]):

D Fe 2+ = 7.2 • 10 -10 m 2 .s -1 , D Cl -= 2 • 10 -9 m 2 .s -1 , D Na + = 1.3 • 10 -9 m 2 .s -1 , D H + = 9.3 • 10 -9 m 2 .s -1 and D Fe(OH) + = 0.8 • 10 -9 m 2 .s -1 .
The initial data used in computation are given by:

x 0 d = 2 mm, C 0 Fe 2+ (x) = 9.8 • 10 -7 mol/L, C 0 H + (x) = C 0 Fe(OH) + (x) = 2.2 • 10 -8 mol/L, C 0 Na + (x) = 10 -3 mol/L and C 0 Cl -(x) ≈ 10 -3 mol/L, for 0 < x < x 0 d ,
where x 0 d is the initial depth of the crevice and C 0 Cl -is given to satisfy the local electroneutrality condition of the solution.

Computation of the minimum value φ m for propagation simulation

The electric potential φ 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 φ m to be imposed depends on the initial data. Let φ min m be such a value. Then, φ min m satisfies the Nernst equation which is given in its dimensional form by:

ε eq := φ min m -φ * = E 0 Fe 2+ /Fe + RT z Fe 2+ F ln C Fe 2+ C ref , C ref = 1 mol/L (4.2)
where ε eq is the electrode potential of iron in contact with the pit/crevice solution at equilibrium, E 0 Fe 2+ /Fe denotes the standard potential for Fe 2+ ↔ Fe, which is known to be E 0 Fe 2+ /Fe = -0.44 V relative to the standard hydrogen electrode at room temperature. The value φ * is usually called φ eq in the literature. It denotes the potential at equilibrium for which the Butler-Volmer flux (2. The region above the red curve corresponds to the range of φ m where dissolution occurs (φ m > φ min m -φ * and f Fe 2+ > 0) and the region below corresponds to the reduction field (φ m < φ min m -φ * ). To be close to equilibrium, the φ m value needs to be close to the red curve (φ m ≈ φ min m -φ * and f Fe 2+ < 0).

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 0 Cl -in the pit/crevice solution.

• the impact of the potential applied to the surface of pure iron φ 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 φ 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 • 10 -9 m 2 /s, L = 10 3 µm = 1 mm, C 0 = 10 -3 mol/L, (4.3)
and for the following physical and chemical parameters [23]:

D Fe 2+ = 7.2 • 10 -10 m 2 .s -1 , D Cl -= 2 • 10 -9 m 2 .s -1 , D Na + = 1.3 • 10 -9 m 2 .s -1 , (4.4
)

k a = 89.1 mol.m -2 .s -1 , k c = 1.2 • 10 -13 mol.m -2 .s -1 (Appendix B of [12, p.195]) (4.5)
with an initial crevice depth x 0 d := x d (0) = 1 mm. 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 4.7 describes the evolution of the potential drop ∆φ = φ(x d (t), t) -φ(0, t) = φ(x d (t), t), which coincides, in this case, with the potential at the bottom of the crevice φ x d (t), t . It shows that the potential drop converges to a certain constant potential value that depends on C 0 NaCl . Physically, this can be described as a convergence to a state of equilibrium after a large time of crevice propagation : C Fe 2+ x d (t), t ≈ constant and φ x d (t), t ≈ constant as t → ∞. 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. 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.

Influence of the iron potential φ m

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

Both graphs illustrate the very high sensitivity of the propagation speed to the metal potential value φ m . Figure 4.8 shows that the evolution of the corrosion speed is constant for the lower values of φ m . Moreover, it is clear that for the higher values of φ m (for example, φ m = -0.2 V), the crevice cannot keep a constant speed of growth.

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:

∂C j ∂t = ∂ ∂x D j ∂C j ∂x + z j F RT C j ∂φ ∂x , t > 0, 0 < x < x d (t), j = 1, . . . , M (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 diff j + J mig j = -D j ∂C j ∂x -D j z j C j ∂φ ∂x , (4.7)
where

J diff j = -D j ∂C j ∂x
is the diffusion flux which is associated with Fick's law and

J mig j = -D j z j C j ∂φ ∂x
is the migration flux caused by the electric field E = -∂φ ∂x . The proportion of transport by migration of the species j versus total transport is computed as follows :

% migration transport j = 100 • |J mig j | |J diff j | + |J mig j | (4.8)
Two cases of localized corrosion were evaluated: crevice corrosion and pit corrosion. For crevice corrosion, an initial crevice depth x 0 d = 100µm and an initial iron concentration C 0 Fe 2+ = 10 -6 mol/L were assumed according to usual crevice configuration. The computations showed that the transport of Cl -and 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 Fe 2+ transport were quite different. Figures 4.12 and 4.13 show the profile of the 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 φ m when the crevice depth reaches 200 µm. These figures show that the Fe 2+ concentration and the potential drop both increase as φ m increases. And, for a fixed value of φ m , lower NaCl concentration leads to lower Fe 2+ concentration (Figure 4.12) and higher potential drop (Figure 4.13). This higher potential drop is no surprise, since lower NaCl and Fe 2+ concentrations imply higher electric resistance of the crevice solution. 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 encoutered in pit configuration but needs to be checked. It is done below.

Parametric study for pit propagation

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 φ 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 mol/L and C Na + (0, 0) = C Cl -(0, 0) = 1 mol/L. ( 4.9)
At the pit bottom (x = x d (0)): In the pit solution (0 < x < x d (0)):

C Fe 2+ (x d (0), 0) = 5 mol/L, C Na + (x d (0), 0) = 1 mol/L, C Cl -(x d (
C Fe 2+ (x, 0) = linear profile from 10 -6 mol/L at x = 0 to 5 mol/L at x = x d (0), C Cl -(x, 0) = linear profile from 1 mol/L at x = 0 to 11 mol/L at x = x d (0), C Na + (x, 0) = 1 mol/L for all 0 < x < x d (0) (4.11)
The repassivation potential φ stable m , 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 φ stable m is shown relative to the NHE reference for several initial pit depths x 0 d := x d (0). As observed, the repassivation potential φ stable m 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]). mass transport restriction results for 1-D pits to sustain the steady-state concentration necessary to corrode stably with a film-free surface. For shallower 1-D pits, the diffusion out of the pit would overwhelm the rate of metal dissolution, leading to repassivation. The view that restricted mass transport impedes repassivation has been substantiated by Sridhar and Cragnolino, 32 as well as Srinivasan et al., 30 who have shown that crevice formation in large-area samples can lead to the measurement of more active repassivation potential values even when the charge density passed is low. Finally, studies using artificial pits 27,28,30 have shown that the depth at which the measured repassivation potential approaches a plateau is approximately eight to ten times the pit diameter, as can be seen in the data reproduced in Figure 15b. These latter results lend credence to the proposition that at this pit depth, a critical steady state concentration (lower than 100% saturation) 30,32,34 is attained at the pit base, with equal rates of both metal dissolution and diffusion out of the pit, thus marking the transition between pit stability and repassivation.

Conclusions

This study employed artificial pit experiments as well as finite element modeling to elucidate the effects of pit depth on the flux emanating from 1-D pits, providing commentary on proper collection, analysis, and interpretation of pit stability data from artificial experiments. For shallow pits, the hemispherical boundary layer at the pit mouth significantly affected the flux characteristics. It was only once pits had attained depths of around eight times the diameter that the overall diffusion length could be approximated well by the pit depth alone. The examination of results from artificial pit experiments reinforced the importance of obtaining data from deep pits so that pit stability parameters extracted not only reflected true 1-D mass transport based on the pit depth, but also were not affected by the bulk electrolyte. Artificial pit experiments performed in several configurations ruled out the influence of convective transport and the precipitation of corrosion products on the mass transport in pits of small diameter. The trends observed in published experimental data were also explained by the inferences of bulk electrolyte effects obtained from this work. Finally, the results of this investigation were also used to rationalize the observed dependence of repassivation potential on pit depth, thereby providing a quantitative connection between pit stability and repassivation studies. 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.

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, φ m is set at 0.1 V/NHE to be higher than the critical φ 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. 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 φ have been removed (but where of course the potential is kept in the Butler-Volmer flux f Fe 2+ (φ, 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+ (φ, C 1 )) at the previous time step t n . Then, we compute at time t n+1 , the profile of the potential φ 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). 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, φ 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.

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. below for the ALE system (A.3)-(A.5). The procedure to determine the initial potential φ 0 i for all i = 1, . . . , 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:

D + ∆x u i = u i+1 -u i ∆x , D - ∆x u i = u i -u i-1 ∆x . (B.5)
• Step 1 : Interface displacement (A.5). We start by computing the new position x n+1 d of the interface using the Butler-Volmer formula

x n+1 d = x n d + v n ∆t with v n = Ω F e f Fe 2+ (φ n N , C 1,n N ) at x = x n d = x n N , (B.6) i.e. with v n x d (t n ) (forward Euler), where f Fe 2+ (φ n N , C 1,n N ) is given by (see (A.7)) f Fe 2+ (φ n N , C 1,n N ) = k a exp φ m -φ n N -k c C 0 C ref C 1,n N exp -( φ m -φ n N ) . (B.7)
choose ∆t small enough so that v n ∆t does not exceed ∆x.

• ] are given by:

x n+1 i := x i = ψ n+1 (x n i ), for i = 0, . . . , N. (B.9)

• Step 3 : Full implicit in time and space discretizations of system (A.3),(A.4): computation of C j,n+1 i and φ n+1 i for j = 1, . . . , M ; i = 0, . . . , N .

We take t = t n+1 in (A.3) at time t = t n+1 . We solve the following problem for C j,n+1 i and φ n+1 i for i = 0, . . . , N :

-The electrodiffusion equations (A.3a) are discretized as : C j (x i , t n ) = C j (ψ n+1 (x n i ), t n ) = C j (x n i , t n ), φ n i φ( x i , t n ) = φ(ψ n+1 (x n i ), t n ) = φ(x n i , t n ), (B.17)

C j,n+1
Hence we see that C j,n i and φ n i are also approximations of C j (x n i , t n ) and φ(x n i , t n ) respectively, where C j and φ are the solutions of the original system (3.1), (2.22), (2.23).

With m = (M + 1) × (N + 1), the nonlinear problem (B.10)-(B.16) forms a nonlinear system of m equations with m unknowns S = {C j,n+1 i } j=1,...,M , φ n+1 i for i = 0, . . . , N . This problem can be written as F(S) = 0 and we use a Newton method to numerically solve it, with S 0 = {C j,n i } j=1,...,M , φ n i ; for i = 0, . . . , 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: 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 = k f k b = C 4 C ref C 5 C ref C 1 C ref C H 2 O C ref 10 -9.316 (C.4)
where C H 2 O = C ref = 1 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: 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 > > 10 5 mol.m -3 .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 mol.m -3 .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].

                 ∂C 1 ∂t = ∂ ∂x D 1 ∂C 1 ∂x + z 1 F RT C 1 ∂φ ∂x - k f C ref C 1 - 1 K 1 C 4 C 5 C ref , (

Figure 2 . 1 :

 21 Figure 2.1: A one dimensional corrosion pit.

j = 1 ,

 1 . . . , M, are the ion fluxes, 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.

4 .

 4 Numerical simulations of the unidimensional pit/crevice 4.1. Comparison of the numerical model with a steady-state model from the literature

Figure 4 .

 4 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 4 .

 4 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 4 . 1 :

 41 Figure 4.1: Concentration profiles along the cavity length for a crevice with passive walls given by Sharland model [6].

Figure 4 . 2 :

 42 Figure 4.2: Evolution of the solution chemistry in the crevice given by the developed model with moving interface for the initial data : φm = 0.04 V |NHE , C 0 Fe 2+ = 10 -6 mol/L, C 0 NaCl = 10 -3 mol/L and x d (0) = 2 mm.

  7) vanishes. The curve of φ min m -φ * as a function of the initial concentration of the iron cations at the bottom of the pit/crevice C 0 Fe 2+ (x 0 d ) is plotted in Figure 4.3 (using equation (4.2)).

Figure 4 .

 4 Figure 4.3: φ min m -φ * as a function of the initial concentration of iron cations at the pit/crevice bottom, given by (4.2).

4. 3 . 1 .

 31 Influence of the initial chloride concentration C 0 Cl -. In this paragraph, simulations are given for a fixed value of φ m = -0.2 V |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 0 NaCl = 1 mol/L is roughly 12 times faster than for C 0 NaCl = 0.001 mol/L.

Figure 4 . 4 :

 44 Figure 4.4: Influence of the initial concentration of chloride ions on the time evolution of the crevice depth.

Figure 4 .

 4 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 ).

  , m = 0.2 V|NHE, C 0 Fe 2 + = 10 6 mol/L, t = 500 hours

Figure 4 . 5 :

 45 Figure 4.5: Evolution of the crevice depth after 500 hours of the crevice propagation for different values of initial chloride concentrations C 0 NaCl (with x 0 = 1 mm, φm = -0.2 V(NHE), C 0 Fe 2+ = 10 -6 mol/L, t = 500 hours).

Figure 4 . 4 ,Figure 4 . 6 :

 4446 Figure 4.6: Influence of the initial concentration of NaCl on the time evolution of the metallic cation concentration C Fe 2+ at the crevice bottom.

Figure 4 .

 4 Figure 4.6 shows that for φ m = -0.2 V| NHE , concentration of Fe 2+ at the bottom of the crevice increases until reaching a quasi-stationary value ≈ 82.2 mol/L. This high concentration

Figure

  Figure 4.7 describes the evolution of the potential drop ∆φ = φ(x d (t), t) -φ(0, t) = φ(x d (t), t), which coincides, in this case, with the potential at the bottom of the crevice φ x d (t), t . It shows that the potential drop converges to a certain constant potential value that depends on C 0 NaCl . Physically, this can be described as a convergence to a state of equilibrium after a large time of crevice propagation : C Fe 2+ x d (t), t ≈ constant and φ x d (t), t ≈ constant as t → ∞. Incidentally, as we see in Figure4.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 4 . 7 :

 47 Figure 4.7: Influence of the initial concentration of NaCl on the time evolution of the potential drop in the crevice solution with time.

Figure 4 . 8 :

 48 Figure 4.8: Influence of the potential applied to the surface of the pure iron φm on the evolution of the crevice depth in time.

Figure 4 .

 4 Figure 4.9 shows the value of the final crevice depth after 500 hours of propagation for different values of the metal potential φ m between -0.35 V| NHE and -0.2 V| NHE .

Figure 4 . 9 :

 49 Figure 4.9: Estimation of the depth evolution law after 500 hours of the crevice propagation as a function of the metal potential φm.

Figure 4 .Figure 4 . 10 :Figure 4 . 11 :

 4410411 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-

Figures 4. 14

 14

and 4 .

 4 15 show the proportion of transport by migration for Fe 2+ species for the same set of values of NaCl concentration and φ m as inFigures 4.12

  and 4.13, when the crevice depth reaches 200 µ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 Fe 2+ transport by migration increases with the metal potential value φ m . Moreover, for a fixed value of φ 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 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 4 . 12 :

 412 Figure 4.12: Profile of Fe 2+ concentration in crevice solution as function of metal potential and NaCl concentration.

Figure 4 . 13 :

 413 Figure 4.13: Potential drop in crevice solution as function of metal potential and NaCl concentration.

Figure 4 . 14 :

 414 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.

Figure 4 . 15 :

 415 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.

  0), 0) = 11 mol/L (by imposed local electroneutrality of the solution). (4.10)

Figure 4 . 16 :Figure 15 .

 41615 Figure 4.16: Identification of the minimal applied potential to ensure the stability of the pit as a function of its initial depth.

Figure 4 . 17 :

 417 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].

Figure 4 . 18 :

 418 Figure 4.18: Profile of Fe 2+ concentration in 100 nm deep pit.Figure 4.19: Potential profile in 100 nm deep pit.

Figure 4 . 19 :

 419 Figure 4.18: Profile of Fe 2+ concentration in 100 nm deep pit.Figure 4.19: Potential profile in 100 nm deep pit.

Figure 4 . 20 :

 420 Figure 4.20: Proportion of migration in transport (%) in 100 nm deep pit.

Figure 4 . 21 :

 421 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 φm = -0.2 V/NHE.

Figure 4 . 22 :

 422 Figure 4.22: Potential profiles in the pit obtained for the pure diffusion model (in yellow) and the diffusion-migration model (in black) for φm = -0.2 V/NHE.

Step 2 :]

 2 Computation of the new mesh of [0, x n+1 d ]. Taking t = t n and t = t n+1 in (A.1) leads to introduce ψ n+1 : [0, x n d ] -→ [0, x n+1 d (N + 1) grid points of the interval [0, x n+1 d

1 --N-Remark 1 .

 11 j = 1, . . . , M and i = 1, . . . , N -Boundary condition at x = x n+1 N = x n+1 d for the metallic cation (see (A.3b)):J 1,n+1 N = -(1 -Ω F e C 1,n+1 N ) f Fe 2+ (φ n+1 N , C 1,n+1N). (B.12) Boundary conditions at x = x n+1 for j = 1, . . . , M .-Boundary conditions at x = x n+1 0 = 0 (see (A.3d), (A.4b)):C j,n+1 0 = β j for j = 1, • • , M such that Electroneutrality condition (see (A.4a)):The convective term in (B.10) is treated with an upwind scheme involving the discrete operator D - ∆x since the velocity flux f F e 2+ is expected to be positive.

Remark 2 .

 2 Due to (B.8), we have that C j,n i

  ∂x + R i (C 1 , . . . , C M ), j = 1, . . . , M. (C.1)where R i (C 1 , . . . , 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:Fe 2+ + H 2 O kf ---kb Fe (OH) + + H + . (C.2)where k f and k b are the forward and backward reaction rate constants. For simplicity, we set C 1 := C F e 2+ , C 2 := C N a + , C 3 := C Cl -, C 4 := C H + , C 5 := C F e(OH) + . (C.3)

  C.5a)

	∂C j ∂t	=	∂ ∂x	D j	∂C j ∂x	+	z j F RT	C j	∂φ ∂x	for j = 2, 3				(C.5b)
	∂C j ∂t	=	∂ ∂x	D j	∂C j ∂x	+	z j F RT	C j	∂φ ∂x	+	k f C ref	C 1 -	1 K 1	C 4 C 5 C ref	for j = 4, 5	(C.5c)
	The equilibrium state of the reaction (C.2) is achieved when k f → +∞, leading to the equilibrium
	relation						C 1 -	1 K 1	C 4 C 5 C ref	= 0.					(C.6)

Acknowledgments

Financial assistance from the Office of the Undersecretary of Defense Technical Corrosion Collaboration via the United States Air Force Academy Contract FA7000-10-2-001 is gratefully acknowledged.

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.

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 0 be fixed. The ALE formulation consists in transforming the current domain [0, x d (t)] for t 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 0:

The function ψ is smooth with respect to t and linear in x. This particular choice of ψ 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

The new functions φ(x, t) and C j (x, t) together with the (unchanged) function x d satisfy the system :

along with the free boundary condition

where

and f Fe 2+ ( φ, C 1 ) is given in (2.17), namely

where the constants φ m , k a and k c in (A.7) are given in (2.18). The initial concentrations C j,0 (x) = C 0 j (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 ∆t > 0 be the discretization time step such that t n = n ∆t be the discrete times for n ∈ N * . The time increment ∆t = t n -t n-1 varies with n. We start with an initial pit bottom position

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

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

given, for i = 0, . . . , N ; j = 1, . . . , M φ 0 0 = 0 and φ 0 N = φ(x 0 d , 0) given.

(B.3)

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

In the ALE scheme, the number (N + 1) of grid points remains constant at each time step t = t n , n 0, but the grid points x n i (i = 0, . . . , N ) of the interval [0, x n d ] 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:


Pierre Barbault

Matthieu Kowalski

Charles Soussen

Estimation du MAP Marginal d'un signal Bernoulli-Gaussien: une approche par relaxation continue

Cet article traite de l'estimation du support de signaux parcimonieux dans le cadre de problèmes inverses. A l'aide d'un a priori Bernoulli-Gaussien, le support est estimé en maximisant la loi a posteriori marginalisée par rapport aux amplitudes des sources. On propose un algorithme d'Espérance-Maximisation (EM), où la variable discrète codant le support est relâchée dans un espace continu. Une évaluation empirique sur signaux simulés avec un opérateur issu de la magnétoencephalographie montre la pertinence de cette approche en comparaison avec l'approche de maximisation de la loi jointe du support et des amplitudes, interprétable comme la minimisation d'un critère avec pénalisation ℓ 0 .

Introduction

On considère le problème inverse linéaire y = Hx + n avec H ∈ R M ×N , x ∈ R N et y ∈ R M . Le vecteur x inconnu est supposé parcimonieux, et modélisé par un processus Bernoulli-Gaussien (BG). Pour tout n, les éléments x[n] sont supposés indépendants, identiquement distribués (i.i.d.), avec

p(x[n]) = p 2πσ 2 x exp - x[n] 2 2σ 2 x + (1 -p)δ(x[n])
où δ est la distribution de Dirac. Le bruit n est supposé indépendant de x, blanc et Gaussien : n ∼ N (0, σ 2 0 I). On écrit x comme produit de deux variables aléatoires indépendantes : ∀n, x[n] = q[n]r[n], où q[n] est une variable binaire représentant le support de x (q[n] = 1 si x[n] ̸ = 0) et r[n] représente les amplitudes du signal. q[n] suit une distribution de Bernoulli : q[n] ∼ B(p) et r[n] une loi gaussienne : r[n] ∼ N (0, σ 2

x ). En notation matricielle, cela donne

x = Qr avec Q = Diag(q). (1)
L'estimateur du MAP de x = {q, r} peut être formulé comme le maximum de la vraisemblance jointe a posteriori de (q, r) ou comme la vraisemblance de la loi marginale a posteriori de q. Ces approches conduisent à deux estimateurs distincts [START_REF] Champagnat | Unsupervised deconvolution of sparse spike trains using stochastic approximation[END_REF] : le MAP Joint (JMAP) et le MAP Marginal (MMAP). L'expression du JMAP est obtenue à partir de p(q, r | y) : . En utilisant le changement de variable (1), l'estimation au sens du JMAP se ramène à un pro-Ce travail est supporté par l'Agence Nationale de la Recherche à travers le projet BMWs (ANR-20-CE45-0018) blème de minimisation ℓ 2 + ℓ 0 [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] :

argmin q,r
x JMAP = argmin x 1 2σ 2 0 ∥y -Hx∥ 2 + 1 2σ 2 x ∥x∥ 2 +ρ∥x∥ 0 . (3)
À q fixé, le minimiseur de (2) selon r admet une forme explicite r(q). En réintroduisant l'expression dans la fonction coût (2), le critère JMAP se récrit :

qJMAP = argmin q∈{0,1} N y T Γ y (q) -1 y + ρ∥q∥ 0 (4) avec Γ y (q) = σ 2 0 I + σ 2 x HQQ T H T . (5)
Une fois le support q JMAP estimé, les amplitudes non nulles sont déduites au sens des moindres carrés :

x(q) = argmax x p(x|y, q) = σ 2 x QH T Γ -1 y (q)y (6)
avec q = q JMAP . L'estimateur du MMAP consiste à minimiser p(q | y), après marginalisation des amplitudes r du signal [START_REF] Champagnat | Unsupervised deconvolution of sparse spike trains using stochastic approximation[END_REF]. C'est un estimateur naturel dans les applications où le support contient l'information principale, comme pour la localisation de sources cérébrales en Magneto/ElectroEncephalographie (M/EEG) [START_REF] Koessler | Source localization of ictal epileptic activity investigated by high resolution EEG and validated by SEEG[END_REF]. Il semble plus approprié d'estimer d'abord le support par le MMAP, qui correspond à l'estimateur de Bayes associé à l'erreur discrète 0 -1, puis de reconstruire les valeurs de r.

Contributions et plan de l'article. On propose en Section 2 une approche par Espérance-Maximisation (EM) dédiée à l'estimation du MMAP. En Section 3, les contraintes d'optimisation binaire sont relâchées dans un domaine continu, [0, 1] N . La solution du problème relâché sert ensuite d'initialisation à chaud de l'algorithme EM pour variables binaires. Enfin, les expériences présentées en Section 4 viennent appuyer la validité de l'approche.

2 Estimation du support via le MMAP L'estimateur du MMAP est basé sur la maximisation de p(q|y) = p(q, r|y) dr sur {0, 1} N . Comme décrit dans [START_REF] Champagnat | Unsupervised deconvolution of sparse spike trains using stochastic approximation[END_REF], il conduit à minimiser

-log p(q|y) = 1 2 y t Γ -1 y (q)y + 1 2 log |Γ y (q)| + ρ∥q∥ 0 . (7)
qui s'avère être un problème NP-difficile. On propose une résolution sous-optimale basée sur un algorithme EM. En suivant [START_REF] Figueiredo | An EM algorithm for waveletbased image restoration[END_REF], le modèle d'observation est décomposé en deux modèles liés à l'opérateur direct et au bruitage additif :

y = Hz + e et z = Qr + b. (8)
On impose les distributions

e ∼ N (0, Γ e ) et b ∼ N (0, σ 2 b I) (9)
définissant des vecteurs indépendants, gaussiens, tels que

Γ e + σ 2 b HH T = σ 2 0 I avec σ 2 b ≤ σ 2 0 ∥HH T ∥ . ( 10
)
où ∥ .∥ est la norme spectrale d'une matrice. Cette décomposition permet de dériver l'algorithme Iterative Shrinkage/Thresholding pour les problèmes de Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]/Basis Pursuit Denoising [START_REF] Chen | A SAUNDERS : Atomic decomposition by basis pursuit[END_REF]. Elle est aussi utilisée dans [START_REF] Kowalski | An unsupervised algorithm for hybrid/morphological signal decomposition[END_REF][START_REF] Barbault | Parameter estimation in sparse inverse problems using Bernoulli-Gaussian prior[END_REF] afin d'estimer les paramètres (ρ, σ 2 x ) d'un modèle Bernoulli-Gaussien. Pour maximiser p(q|y), considérons z comme variable cachée admissible. L'algorithme EM s'écrit :

q t+1 = argmin q∈{0,1} N Q(q, q t ) Q(q, q t ) = E z|y,q t[-log p(z, q|y)] . (11)
L'étape E et l'étape M sont développées dans la suite.

Etape-E

Dans un but de concision, la notation suivante sera adoptée :

f (q) + κ κ = f (q)
pour faire référence à toute égalité à une constante κ près (qui ne dépend pas de q). Par la linéarité de l'espérance et le fait que p(q|z, y) = p(q|z), on a

Q(q, q t ) κ = E z|y,q t[-log p(q|z)] (12)
où κ = E z|y,q t[-log p(z|y)]. La règle de Bayes mène à :

Q(q, q t ) κ = E z|y,q t[-log p(z|q)] + ρ∥q∥ 0 (13) avec ρ défini en (2), avec maintenant κ = E z|y,q t[log p(z) -log p(z|y)]. En utilisant Eqs. (8) à (10), on a directement z|q ∼ N (0, Γ z (q)) avec Γ z (q) = σ 2 b I + σ 2 x Q T Q. ( 14
)
La règle de Bayes donne la loi a posteriori de z z|y, q ∼ N (ẑ, Σ)

dont les paramètres sont

ẑ = Γ z (q)H T Γ -1 y (q)y Σ = Γ z (q) -Γ z (q)H T Γ -1 y (q)HΓ z (q). (16)
C'est pourquoi l'espérance s'écrit

E z|y,q t[-log p(z|q)] κ = 1 2 E z|y,q t z T Γ -1 z (q)z + 1 2 log |Γ z (q)| κ = 1 2 (ẑ t ) T Γ -1 z (q)(ẑ t ) + 1 2 Trace[Γ -1 z (q)Σ t ] + 1 2 log |Γ z (q)|.
Enfin, l'étape-E s'écrit comme la somme séparable

Q(q, q t ) κ = 1 2 N n=1 Q n (q[n]) (17)
avec

Q n (q[n]) = (ẑ t [n]) 2 + Σ t [n, n] σ 2 b + σ 2 x q[n] 2 +log(σ 2 b +σ 2 x q[n] 2 )+2ρq[n]
Ainsi, la minimisation de Q se réduit à N problèmes de minimisation 1D, i.e., la minimisation des critères Q n .

Etape-M

Introduisons la notation

γ t n = (ẑ t [n]) 2 + Σ t [n, n]. La minimi- sation de Q n s'écrit : q t+1 [n] = argmin q∈{0,1} γ t n σ 2 b + σ 2 x q 2 + log(σ 2 b + σ 2 x q 2 ) + 2ρq (18)
Après calculs, la solution est q t+1 [n] = 1 si et seulement si

γ t n ≥ σ 2 b σ 2 x + σ 2 b σ 2 x log σ 2 b + σ 2 x σ 2 b + 2ρ . ( 19
)
Cette opération n'est pas sans rappeler la formule de seuillage dans le cadre du problème de débruitage en MAP-Marginal (avec γ t n en lieu et place de z[n] ci-dessous).

Proposition 2.1 ([13, 9, 1]). Soit z = Qr + b, avec les lois a priori définies plus haut pour (q, r). L'estimateur du MMAP de q|z est donné explicitement, pour tout n, par

q[n] = 1 ssi P (q[n] = 1|z[n]) ≥ P (q[n] = 0|z[n]) c.à.d. ssi z[n] 2 ≥ 2σ 2 b σ 2 x + σ 2 b σ 2 x ρ + 1 2 log σ 2 b + σ 2 x σ 2 b .

Résumé de l'algorithme

La méthode complète est résumé dans Alg. 1. Lorsque H = I, e peut être mis à 0 dans (8), et donc σ 2 b = σ 2 0 . Dans ce cas, les matrices de covariance Γ z et Γ y sont égales, et l'algorithme se réduit au cas de l'estimateur du MMAP en débruitage (avec z = y dans Prop. 2.1). Dans la section suivante, une relaxation continue du problème MMAP est proposée, où les variables binaires q[n] sont remplacées par des variables dans [0, 1]. L'idée est d'obtenir une solution « à chaud » qui permettra de mieux initialiser Alg. 1.

Algorithme 1 : Estimateur MMAP du support q.

Entrées : t = 0, q (t) ∈ [0, 1] N , ρ > 0; Sortie : q ∈ {0, 1} N ; while not converged do

Q (t) = Diag(q); Γ (t) z = σ 2 b I + σ 2 x Q (t) Q (t) T ; Γ (t) y = σ 2 0 I + σ 2 x HQ (t) Q (t) T H T ; Σ (t) = Γ (t) z -Γ (t) z H T Γ (t) -1 y HΓ (t) z ; ẑ(t) = Γ (t) z H T Γ (t) -1 y y; for n = 1, . . . , N do γ (t) n = ẑt [n] 2 + Σ t [n, n]; if γ (t) n ≥ 2σ 2 b σ 2 b +σ 2 x σ 2 x ρ + log 1 + σ 2 x σ 2 b then q (t+1) [n] = 1; else q (t+1) [n] = 0; end end t = t + 1 end

Optimisation par relaxation continue

On relâche le problème d'optimisation binaire [START_REF] Lim | Sparse EEG/MEG source estimation via a group lasso[END_REF] en un problème d'optimisation continue. Le domaine de définition {0, 1} N de la fonction coût est remplacé par [0, 1] N . Comme la fonction coût Q(q, q t ) est la somme découplée [START_REF] Wipf | Sparse Bayesian learning for basis selection[END_REF], l'étape M reste un problème 1D semblable à (18).

On remarque que pour q[n] ∈ {0, 1}, alors q[n] 2 = q[n]. Ainsi, minimiser Q n sur {0, 1} revient à minimiser :

γ (t) n σ 2 b + σ 2 x q[n] 2 + log(σ 2 b + σ 2 x q[n] 2 ) + 2ρq[n] 2 .
Considérons la relaxation continue de q[n] sur [0, 1] :

q t+1 n = argmin q∈[0,1] γ t n σ 2 b + σ 2 x q 2 + log(σ 2 b + σ 2 x q 2 ) + 2ρq 2 . ( 20
)
Il s'agit d'un problème 1D avec des contraintes de borne. Le minimiseur non contraint peut être obtenu en cherchant les racines de la dérivée première par rapport à u = q 2 , qui s'écrit

-σ 2 x γ t n (σ 2 b + σ 2 x u 2 ) 2 + σ 2 x σ 2 b + σ 2 x u + 2ρ = 0 . (21)
Lorsque ρ > 0, il existe une seule racine positive :

u[n] = 1 4ρ 1 + 8ρ γ t n σ 2 x -1 - σ 2 b σ 2 x . (22)
Une étude minutieuse de la fonction coût dans (20) (considérée comme une fonction de u = q 2 ) révèle que la dérivée seconde est positive pour le point u[n]. En prenant en compte que q ∈ [0, 1], la solution devient :

q t+1 [n] =      0 si u[n] ≤ 0, 1 si u[n] ≥ 1, u[n] sinon. ( 23
)
Algorithme 2 : Estimation of continuous relaxed q Entrées : t = 0, q (t) ∈ [0, 1] N ; Sortie : q ∈ [0, 1] N ; while not converged do

Q (t) = Diag(q) ; Γ (t) z = σ 2 b I + σ 2 x Q (t) Q (t) T ; Γ (t) y = σ 2 0 I + σ 2 x HQ (t) Q (t) T H T ; Σ (t) = Γ (t) z -Γ (t) z H T Γ (t) -1 y HΓ (t) z ; ẑ(t) = Γ (t) z H T Γ (t) -1 y y; for n = 1, . . . , N do γ (t) n = ẑt [n] 2 + Σ t [n, n]; u[n] = 1 4ρ 1 + 8ρ γ t n σ 2 x -1 - σ 2 b σ 2 x ; q[n] = max{min{ u[n], 1}, 0} ; end t = t + 1; end L'
algorithme obtenu Alg. 2 est proche de l'EM Sparse Bayesian Learning (SBL) de [START_REF] Wipf | Sparse Bayesian learning for basis selection[END_REF], exception faite de l'estimation de q. Comme x|q ∼ N (0, σ 2

x Q), lorsque q est continu, on retrouve la loi a priori utilisée dans SBL. La relaxation proposée peut être vue comme une approche SBL où les variances sont bornées par [0, σ 2 x ] N . Enfin, quand ρ = 0, on a p = 1 2 dans le modèle BG. Les mises à jour (20) de q[n] sont données par

q[n] 2 = min 1, γ t n -σ 2 b σ 2 x + (24)
qui conduit exactement à l'EM SBL [START_REF] Wipf | Sparse Bayesian learning for basis selection[END_REF] si σ 2 x = 1 et σ 2 b = 0. Concrètement, Alg. 2 est lancé pour différentes valeurs de l'hyperparamètre ρ. Les résultats sont utilisés comme initialisation à chaud de l'EM binaire pour le MMAP (Alg. 1) avec les vraies valeurs des paramètres du modèle. Cette initialisation permet d'éviter que Alg. 1 converge trop rapidement vers un mauvais minimum local.

Expériences

Le signal Bernoulli-Gaussien (q, r) est simulé avec les paramètres p = 0.05 et σ 2

x = 1, et le bruit blanc gaussien est de variance σ 2 0 = 0.01 ( ce qui correspond à un SNR d'entrée théorique d'environ 7 dB ). L'opérateur H utilisé est une sous-matrice issue de l'opérateur direct (leadfield) en magnéto-encéphalographie (M/EEG), de taille 272 × 600. Cette matrice contient des colonnes hautement corrélées. Elle est représentative de la difficulté des problèmes inverses en M/EEG. L'estimation de q au sens du MMAP est comparée aux estimations correspondantes x(q) telles que dans [START_REF] Figueiredo | An EM algorithm for waveletbased image restoration[END_REF], ainsi qu'au classique MAP Joint de q et x obtenu en minimisant (3). Pour ces méthodes, on utilise un algorithme de descente proximale similaire à l'algorithme d'Iterative Hard Thresholding (IHT) [START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF] avec redémarrage à chaud. Les comparaisons sont faites selon les métriques suivantes : le taux de Faux Positifs (erreur de Type-I), le taux de Faux Négatifs (erreur de Type-II), le SNR des x(q) estimé selon [START_REF] Figueiredo | An EM algorithm for waveletbased image restoration[END_REF]

-3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3 -3.3
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1. Les erreurs de Type-I et de Type-II sont représentées en fonction de la taille du support prédit q en Fig. 1. On observe que la solution du MMAP est de meilleure qualité que celle du JMAP, que ce soit en terme d'erreur de type-I ou II et de SNR (surtout, pour les métriques Type-II et RSB). De plus, l'algorithme se comporte bien pour des opérateurs H très corrélés. D'autres expérimentations (non montrées ici) montrent que pour des matrices H gaussiennes, modérément corrélées, les estimateurs du JMAP et du MMAP donnent des résultats similaires. La version Marginale reste plus gourmande en temps de calcul à cause de la présence de l'inverse de la matrice Σ t dans Algs. 1 et 2.

Discussion et conclusion

Cet article propose une approche par MMAP pour sélectionner le meilleur support possible, car il s'agit de l'estimateur bayésien pour l'erreur 0 -1. Ce problème de sélection de support peut être vu commme un problème de sélection de variables très étudié en statistique [START_REF] O'hara | SILLANPÄÄ : A review of Bayesian variable selection methods : What, how and which[END_REF], aussi connu comme best subset selection problem [START_REF] Hastie | TIBSHIRANI : Best subset, forward stepwise or lasso ? analysis and recommendations based on extensive comparisons[END_REF]. En classification, ce problème est en général étudié dans un contexte de « prédiction » de la classe de nouvelles données ; ainsi, l'estimation de x est primordiale. Dans [START_REF] Hastie | TIBSHIRANI : Best subset, forward stepwise or lasso ? analysis and recommendations based on extensive comparisons[END_REF], une comparaison exhaustive est faite entre le Lasso et les solutions optimales ℓ 0 en utilisant un solveur MILP venant de [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming : Formulations and computational performance[END_REF]. Dans ce contexte, les résultats des approches ℓ 0 et Lasso sont similaires. Cependant, la matrice utilisée dans ces simulations n'est pas aussi mal posée que l'opérateur de M/EEG. Notre étude permet d'entamer une réflexion sur les relations existantes entre les minimiseurs du MAP Joint et du MMAP. Nos futurs travaux s'attaqueront à l'impact de la relaxation continue. De plus, une comparaison en profondeur de ces méthodes avec les approches MCMC appliquées aux a priori BG [START_REF] Ge | CARPENTIER : Enhanced sampling schemes for mcmc based blind Bernoulli-Gaussian deconvolution[END_REF] et les EMVS proposés par [START_REF] Rovcková | EMVS : The EM approach to Bayesian variable selection[END_REF] (basés sur des mélanges de gaussiennes, et donc s'approchant de z dans (8)). Du point de vu pratique, la méthode sera particularisée au cas de la structure spatio-temporelle des signaux issus de M/EEG [START_REF] Lim | Sparse EEG/MEG source estimation via a group lasso[END_REF].

2 x ∥r∥ 2 +

 22 ρ∥q∥ 0 , (2) avec ρ = log 1-p p

2 ||x

 2 et calculé en prenant SN R = 10 log 10 ||x| 2

FIGURE 1 :

 1 FIGURE 1 : Estimateurs du MAP marginal vs du MAP joint, comparés pour p = 0.01 et σ 2 0 = 0.01. La mesure de SNR des signaux reconstruits x est donnée pour certains points (de même que le score du meilleur SNR atteint).