# Thermodynamically consistent modelling of the corrosion of iron in the context of deep subsurface nuclear waste repositories

Clément Cancès, Claire Chainais-Hillairet, Benoît Merlet, Federica Raimondi, Juliette Venel

► To cite this version:

Clément Cancès, Claire Chainais-Hillairet, Benoît Merlet, Federica Raimondi, Juliette Venel. Thermodynamically consistent modelling of the corrosion of iron in the context of deep subsurface nuclear waste repositories. EURAD joint program. 2023. hal-04210782

HAL Id: hal-04210782

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

Submitted on 19 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

European Joint Programme
on Radioactive Waste Management

Thermodynamically consistent modelling of the
corrosion of iron in the context of deep subsurface
nuclear waste repositories

¹Cancès, C., ¹Chainais-Hillairet, C., ¹Merlet, B., ²Raimondi, F. and ³Venel, J.

¹Univ. Lille, CNRS, Inria, UMR 8524 - Laboratoire Paul Painlevé, F-59000 Lille, France

²Department of Mathematics, Università di Salerno, 84084, Fisciano, SA, Italy

³Univ. Polytechnique Hauts-de-France, CERAMATHS, FR CNRS 2037,
F-59313 Valenciennes, France

The project leading to this application has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 847593.

<figure><img src="image_3.png" /></figure>
Abstract

Corrosion occurring in subsurface nuclear waste repositories yields the production of dihydrogen (H₂), the possible accumulation of which being a major security concern. The development of accurate and theoretically assessed mathematical and numerical models is therefore a priority to quantify the *in situ* production of H₂ along large time scales. Despite important efforts of the mathematical community during the last 15 years, there is so far no satisfactory mathematical framework for the so-called Diffusion Poisson Coupled Model (DPCM) proposed by Bataillon and collaborators (Bataillon et al., 2012, 2010). This model describes the evolution of the oxide layer covering the metal by taking into account the oxidation of the metal, the transfer within the oxide layer of the charge carriers driven by some self-consistent electric potential, and the dissolution of the oxide at the interface with some aqueous solution.

The main reason for the aforementioned gap in the theory is that no thermodynamic potential, serving as a Lyapunov functional, has been shown to be dissipated along time, in accordance to the second principle of thermodynamics. Assuming pressure and temperature to be constant, Gibbs free energy is indeed expected to decay, up to some exchange with the surrounding metal and solution. We propose an update of the DPCM model which fulfils some variant of Onsager's reciprocal relation, ensuring therefore the compatibility of the model with the second principle of thermodynamics.

The main differences with the original DPCM model are the following: (i) the transport of ferric cations and oxygen vacancies in the oxide are driven by a vacancy diffusion process with a nonlinear mobility including a saturation effect; (ii) the motion of the interfaces is driven by the difference of grand potential (or Landau free energy) density rather than by a difference of chemical potential as in the original DPCM; (iii) a correction in the charge carrier fluxes is incorporated to take the volume expansion stemming from the iron oxidation into account; (iv) an update of the boundary condition for the electrons between the oxide and the metal.

<figure><img src="image_4.png" /></figure>
**Significance Statement**

Corrosion occurring in subsurface nuclear waste repositories yields the production of dihydrogen ($H_2$), the possible accumulation of which being a major security concern. The development of accurate and theoretically assessed mathematical and numerical models is therefore a priority to quantify the *in situ* production of $H_2$ along large time scales. Despite important efforts of the mathematical community during the last 15 years, there is so far no satisfactory mathematical framework for the so-called Diffusion Poisson Coupled Model (DPCM) proposed by Bataillon and collaborators (Bataillon et al., 2012, 2010). This model describes the evolution of the oxide layer covering the metal by taking into account the oxidation of the metal, the transfer within the oxide layer of the charge carriers driven by some self-consistent electric potential, and the dissolution of the oxide at the interface with some aqueous solution.

The main reason for the aforementioned gap in the theory is that no thermodynamic potential, serving as a Lyapunov functional, has been shown to be dissipated along time, in accordance to the second principle of thermodynamics. Assuming pressure and temperature to be constant, Gibbs free energy is indeed expected to decay, up to some exchange with the surrounding metal and solution. We propose an update of the DPCM model which fulfils some variant of Onsager's reciprocal relation, ensuring therefore the compatibility of the model with the second principle of thermodynamics.

The advantage of having some stability encoded by the second principle of thermodynamics is manyfold. First, this stability property is essential to establish the well-posedness of the system of partial differential equations (PDEs) derived along the modelling process for arbitrary large time. It also paves the way to the design of stable numerical methods that remain accurate even for large time horizons.

<figure><img src="image_5.png" /></figure>
Table of contents

<table><tr><td>Abstract</td><td>2</td></tr><tr><td>Significance Statement</td><td>3</td></tr><tr><td>List of figures</td><td>5</td></tr><tr><td>1. Modelling corrosion for nuclear waste repository safety</td><td>6</td></tr><tr><td>1.1 General context and motivation</td><td>6</td></tr><tr><td>1.2 State of the art and scientific positioning</td><td>6</td></tr><tr><td>2. A thermodynamically consistent DPCM model</td><td>7</td></tr><tr><td>2.1 Model description</td><td>7</td></tr><tr><td>2.2 Onsager's reciprocal relation and free energy dissipation</td><td>10</td></tr><tr><td>3. A reduced model with fixed interfaces</td><td>11</td></tr><tr><td>3.1 Model simplification</td><td>11</td></tr><tr><td>3.2 Mathematical analysis</td><td>11</td></tr><tr><td>3.3 Numerical approximation</td><td>12</td></tr><tr><td>4. Conclusion and prospects</td><td>13</td></tr><tr><td>4.1 Conclusion</td><td>13</td></tr><tr><td>4.2 Prospects</td><td>14</td></tr><tr><td>Code source</td><td>14</td></tr><tr><td>References</td><td>15</td></tr></table>

<figure><img src="image_6.png" /></figure>
List of figures

Figure 1 – Schematic representation of the oxide layer................................................................................7

Figure 2 – Schematic representation of the self-consistent electric potential in the oxide ........................8

Figure 3 – Comparison of the current-voltage characteristics corresponding to the new thermodynamically consistent reduced model (vDPCM, orange solide line) and to the original reduced model (DPCM, dashed blue line) for pH = 7 (left), pH = 8,5 (center) and pH = 10 (right). ........................12

Figure 4 – Profiles of the scaled densities of cations (top), electrons (middle) and the electrostatic potential (bottom, in Volts, the physical unit) at the steady state for two values of the applied potential V, that are -0.4 Volts (left) and 0.3 Volts (right). ..................................................................................13

<figure><img src="image_7.png" /></figure>
# 1. Modelling corrosion for nuclear waste repository safety

## 1.1 General context and motivation

The strategy of the French agency for radioactive waste management (ANDRA) for the long term storage of long life and high activity nuclear wastes relies on the subsurface storage in deep geological layers. Galleries dug in a clay geological layer should host iron canisters containing vitrified wastes. Assessing the security of the whole system on large time periods – up to millennials – is a challenge involving a large scientific community gathering researchers and engineers.

In such a scenario, the corrosion of the iron canisters will lead to the generation of dihydrogen ($H_2$), the accumulation of which being a source of explosion risk. The precise quantification of the $H_2$ source over long time period is therefore a scientific priority. It is the motivation to the work reported here, that was carried out in the framework of our research project supported by the EURAD program.

## 1.2 State of the art and scientific positioning

There has been an important effort by ANDRA to develop a model and a numerical tool to simulate in an accurate way the evolution of the magnetite layer protecting the iron canister to be stored in geological repositories. The code CALIPSO, the main developer of which being Christian Bataillon, is currently the reference code for the simulation of the complex model – referred to in the literature as DPCM for diffusion coupled Poisson model – described in (Bataillon et al., 2012, 2010).

The continuous model addressed here consists in one-dimensional convection diffusion equations for the charge carriers coupled with a self-consistent electric potential governed by a Poisson equation. Non-trivial boundary conditions across the interfaces strengthen the coupling. Moreover, the geometry of the oxide layer evolves along time due to the oxidation of the model and to the dissolution of the oxide (see section 2.1 for a similar construction).

From a numerical point of view, the method (Bataillon et al., 2012) implemented in CALIPSO relies on finite volumes for the space discretization, and on a fully coupled backward Euler scheme for the time discretization. A similar strategy is adopted in our work.

Even though CALIPSO gives good results and looks to be stable, only very few elements of mathematical analysis have been provided so far. Let us mention for instance (Chainais-Hillairet and Lacroix-Violet, 2014, 2012) where solutions to a reduced model with fixed interface are studied in the steady and transient case respectively. In the case of a mobile domain, the existence of pseudo-stationary (or travelling-wave) solution has been studied in (Chainais-Hillairet and Gallouët, 2016) under an electric neutrality assumption for the oxide layer, and in (Breden et al., 2021) thanks to a computer assisted proof, but for a single set of parameters. The lack in mathematical foundations for the CALIPSO code stems from the lack of theoretically proven stability for the DPCM model.

Whereas the stability of scalar partial differential equations can take many forms, it becomes much more complex when it comes to highly coupled systems. In such a case where elementary calculations can no longer be carried out, it is now well understood since a few decades that the stability should come from the second principle of thermodynamics. Indeed, thermodynamical potentials like the Gibbs free energy provide natural candidates as Lyapunov functionals. However, proving such a decay requires a very specific structure for the model, which is in general not fulfilled by models built by the apposition of physical laws to model distinct phenomena. A global approach to capture off-diagonal effects is needed. Furthermore, the boundary conditions and the bulk equations should be set accordingly. The DPCM model was not derived in this spirit. Our goal in this project was to update it so that it becomes compatible with thermodynamics, in the sense that the free energy decays along time.
## 2. A thermodynamically consistent DPCM model

### 2.1 Model description

In the same spirit as (Bataillon et al., 2012, 2010), we are interested in a one-dimensional model, the space variable x being the dimensionless rescaled<fnref n="1" /> width of the oxide layer made of magnetite (Fe₃O₄) protecting the metal. In our setting, magnetite is thought of as a one-dimensional lattice made of one fixed ferric cation located on a tetrahedral site surrounded by 4 oxygen anions which are allowed to jump from one octahedral site to a neighbouring one. Two additional mobile ferric cations Fe³⁺ occupying octahedral sites may jump to available neighbouring sites following a vacancy diffusion process too. We also consider the displacement of electrons e⁻ in the conduction band. Rather than the number of oxygen anions O²⁻ per elementary oxide pattern, we keep track of the concentration of oxygen vacancies V²⁺, the number of which per elementary pattern being equal to 4 minus the number of oxygen anions. We denote by u₁, u₂ and u₃ the concentration of ferric cations, electrons and oxygen vacancies respectively, while z₁ = +3, z₂ = -1 and z₃ = +2 denote their respective (rescaled) charge. We also denote by $\bar{u}_1$ = 2 and $\bar{u}_3$ = 4 the maximal number of mobile ferric cations and oxygen vacancies per elementary magnetite pattern.

<figure><img src="image_8.png" /><figcaption>Figure 1 – Schematic representation of the oxide layer</figcaption></figure>

The position of the interface between the oxide and the solution (resp. the metal) is denoted by X₀ (resp. X₁) and is assumed to evolve along time. Therefore, X₀ and X₁ are unknowns of our problem, as well as the self-consistent electrostatic potential Ψ which is deduced from the charge carrier concentrations thanks to a Poisson equation (Bataillon et al., 2010).

More precisely, we set

$$ -\lambda^2 \partial_{xx} \Psi = \rho_{hl} + \sum_{i=1}^{3} z_i u_i \quad \text{in } (X^0, X^1) \qquad (1) $$

where $\rho_{hl} = -5$ stands for the charge of the host lattice made of one iron atom at the tetrahedral site together with 4 oxygens at octahedral sites. The rescaled Debye length is denoted by $\lambda$. The above Poisson equation is complemented with Robin boundary conditions modelling the interfaces as capacitors. More precisely, we set

$$ -\gamma^0 \partial_x \Psi(t, X^0(t)) + \Psi(t, X^0(t)) - \Delta \Psi_{\text{pzc}}^0 = 0 \qquad (2) $$

and

$$ \gamma^1 \partial_x \Psi(t, X^1(t)) + \Psi(t, X^1(t)) - V + \Delta \Psi_{\text{pzc}}^1 = 0. \qquad (3) $$

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

<footnote n="1">¹ The nondimensionalization of the system is not discussed in this contribution. We refer to (Breden et al., 2021) for details on that purpose.</footnote>
In the above formula, $\gamma^0$ and $\gamma^1$ are positive parameters involving the capacity of the interfaces and the rescaled Debye length, whereas V is the electric potential in the metal (the reference potential in the solution is set to 0), and $\Delta\Psi_{pzc}^0$ (resp. $\Delta\Psi_{pzc}^1$) is the voltage drop due to the charge accumulation in the Helmholtz layer at the solution/oxide (resp. oxide/metal) interface, see Figure 2. In what follows, we make use of the shorten notation $\Psi^0$ and $\Psi^1$ for $\Psi(t, X^0(t))$ and $\Psi(t, X^1(t))$ respectively.

<figure><img src="image_10.png" /><figcaption>Figure 2 – Schematic representation of the self-consistent electric potential in the oxide</figcaption></figure>

The equations governing the evolution of the electric potential are fully similar to what was proposed in (Bataillon et al., 2012, 2010). This is also the case of the equations governing the flux of cations and oxygen vacancies across the interfaces, which are prescribed by Butler-Volmer laws. More precisely, at the solution/oxide interface, we set

$$F_i^0 = -k_i^0 u_i^0 e^{z_i \beta_i^0 \Psi^0} m_i^0 + (\bar{u}_i - u_i^0) e^{-z_i \alpha_i^0 \Psi^0}, \quad i = 1 \text{ and } i = 3, \quad (4)$$

with $k_i^0$ and $m_i^0$ being positive parameters, and $\alpha_i^0$ and $\beta_i^0$ being nonnegative parameters summing to 1. Similarly, for the flux of iron atoms through the interface $X^1$, we have

$$F_1^1 = m_1^1 u_1^1 e^{z_1 \beta_1^1 (\Psi^1 - \nu)} k_1^1 - (\bar{u}_1 - u_1^1) e^{-z_1 \alpha_1^1 (\Psi^1 - \nu)}, \quad (5)$$

There is no flux of oxygen atoms through the oxide/metal interface. Translating this fact in terms of flux of oxygen vacancies yields the following formula linking the vacancy flux with the interface motion:

$$F_3^1 = -\bar{u}_3 \dot{X}^1. \quad (6)$$

Denoting by $a_2^P(t) = a_2(u_2(t, X^P(t)))$ the activity of the electrons at the interface $X^P$, $p \in \{0,1\}$, then the electronic exchange with the metal writes

$$F_2^1 = m_2^1 a_2^1 e^{z_2 \beta_2^1 (\Psi^1 - \nu)} - k_2^1 a_2^{\text{met}} e^{-z_2 \alpha_2^1 (\Psi^1 - \nu)}, \quad (7)$$

with $a_2^{\text{met}}$ standing for the activity of the electrons in the metal. The aforementioned expression slightly differs from the one proposed in (Bataillon et al., 2010), which however suggests to set the parameters $\beta_2^1 = 0$ and $\alpha_2^1 = 1$.

Still in (Bataillon et al., 2010), it is suggested that the flux of electrons between the aqueous domain and the oxide has two origins. The ferrous release (FR), which accounts for the reduction of ferric into ferrous cations in the solution, is described by the Butler-Volmer type relation

$$F_2^{0,\text{FR}} = m_2^{0,\text{FR}} e^{-z_2 \alpha_2^{0,\text{FR}} \Psi^0} - k_2^{0,\text{FR}} a_2^0 e^{z_2 \beta_2^{0,\text{FR}} \Psi^0}, \quad (8)$$

for positive $m_2^{0,\text{FR}}$ and $k_2^{0,\text{FR}}$ respectively depending on the activity of the ferrous and ferric cations in the solution. Besides, proton reduction (PR) also induces a flux of electrons governed by

$$F_2^{0,\text{PR}} = m_2^{0,\text{PR}} e^{-z_2 \alpha_2^{0,\text{PR}} \Psi^0} - k_2^{0,\text{PR}} a_2^0 e^{z_2 \beta_2^{0,\text{PR}} \Psi^0}, \quad (9)$$

with $m_2^{0,\text{PR}}$ and $k_2^{0,\text{PR}}$ respectively depending on the activity of dissolved $H_2$ in the solution and on the pH. The electron flux across $X^0$ is then given by

$$F_2^0 = F_2^{0,\text{FR}} + F_2^{0,\text{PR}}. \quad (10)$$

<figure><img src="image_11.png" /></figure>
Denote by R<sub>PB</sub> the Pilling-Bedworth ratio, i.e. the ratio between one elementary cell of oxide over the volume occupied by the required number of iron atoms to form such a pattern (3 for the magnetite). In our case, this ratio is close to 2 for the magnetite and in particular greater than 1, meaning that the oxide has a smaller concentration of iron atoms than the metal. The inertial frame of the oxide then moves with velocity

$$v_{ox} = (1 - R_{PB}) \dot{X}^1. \qquad (11)$$

The conservation of the species i ∈ {1,2,3} then writes

$$\partial_t u_i + \partial_x J_i = 0, \qquad (12)$$

where J<sub>i</sub> denotes the flux of species i in the inertial frame of reference. Owing to (11), it relates to the flux G<sub>i</sub> in the inertial frame of the oxide through formula

$$J_i = G_i + u_i v_{ox}. \qquad (13)$$

The ferric cations and oxygen vacancies bulk fluxes G<sub>i</sub>, i ∈ {1,3}, have nonlinear convection and linear diffusion

$$G_i = -d_i (\partial_x u_i + \eta_i(u_i) z_i \partial_x \Psi), \qquad \eta_i(u_i) = \frac{u_i(\bar{u}_i - u_i)}{\bar{u}_i}, \quad i \in \{1,3\}, \qquad (14)$$

where d<sub>i</sub> > 0 stands for the diffusion coefficient of the species i. Define by

$$\mu_i = \mu_i^* + \log \frac{u_i}{\bar{u}_i - u_i}, \quad i = 1 \text{ and } i = 3, \qquad (15)$$

the chemical potentials of the ferric cations and of the oxygen vacancies – the reference chemical potential μ<sub>i</sub><sup>*</sup> depends on the temperature and on the pressure only, which are assumed to be given constants throughout this report –, and by

$$\xi_i = \mu_i + z_i \Psi \qquad (16)$$

their electrochemical potential. Then the ferric cation and oxygen vacancy fluxes G<sub>i</sub> rewrite

$$G_i = -d_i \eta_i(u_i) \partial_x \xi_i, \quad i \in \{1,3\}. \qquad (17)$$

The nonlinear convection in (14) is characteristic of a vacancy diffusion process. In contrast, there is no limit on the number of electrons in the band conduction provided enough energy is afforded. Denote by a<sub>2</sub> = e<sup>μ<sub>2</sub> - μ<sub>2</sub><sup>*</sup></sup> the activity of the electrons, which is supposed to depend on u<sub>2</sub> only but neither on u<sub>1</sub> nor u<sub>3</sub>, then the flux of electrons in the frame of the oxide is given by

$$G_2 = -d_2 (\partial_x a_2 + z_2 a_2 \partial_x \Psi) = -d_2 a_2 \partial_x \xi_2. \qquad (18)$$

As pointed out in (Bataillon et al., 2012), the boundary fluxes F<sub>i</sub><sup>p</sup> relate to the bulk fluxes thanks to the relation

$$F_i^p = J_i^p - u_i^p \dot{X}^p, \quad i \in \{1,2,3\}, \ p \in \{0,1\}. \qquad (19)$$

In the case of interest where a Boltzmann statistics is used for the electrons, i.e. μ<sub>2</sub> = log(u<sub>2</sub>), then one recovers the same linear drift-diffusion flux (18) as in (Bataillon et al., 2010) up to the term corresponding to the change of inertial frame in (13). This is not the case for the ferric cations and the oxygen vacancies as, additionally to the correction corresponding to the change of frame, convection in our model is nonlinear (14), in opposition to the linear convection postulated in (Bataillon et al., 2012, 2010). These corrections to the reference model (Bataillon et al., 2012) of the state of the art already yield important differences that are highlighted in (Cancès et al., 2023a) when the evolution of the oxygen vacancies – and thus the motion of the interfaces – is neglected.

We introduced a more conceptual difference to handle the evolution of the geometry of the oxide layer and the motion of the interfaces X<sup>0</sup>, X<sup>1</sup>. While the displacement of the interfaces was implicitly driven by
the difference of some electrochemical potential in (Bataillon et al., 2010), this is no longer the case in our model. For thermodynamic reasons that will be made explicit in the next section, we postulate that the displacement of the interface is driven by the jump in the density of the grand potential (or Landau free energy). Denote by $A_i(u_i)$ the chemical contribution to the free energy of species i, so that $dA_i = \mu_i du_i$, then the grand potential is defined by

$$ \Pi = \sum_{i=1}^{3} u_i \mu_i - A_i(u_i) - \rho_{hl} \Psi - \frac{\lambda^2}{2} |\partial_x \Psi|^2, \quad (20) $$

then one readily checks that $\partial_x \Pi = \sum_{i=1}^{3} u_i \partial_x \xi_i$. In particular, $\Pi$ depends on all the species and not only on the oxygen vacancies density, in opposition to the original model (Bataillon et al., 2010). In our model, we suppose that

$$ \dot{X}^1 = \kappa^1 Z^1 (\Pi(X^1) - \bar{u}_3(\mu_3(X^1) + z_3 \Psi(X^1)) - \Pi^{\text{met}}), \quad (21) $$

$$ \dot{X}^0 = v_d^0 + (1 - R_{PB}) \dot{X}^1. \quad (22) $$

The dissolution velocity of the oxide $v_d^0$ is assumed to have the form

$$ v_d^0 = -\kappa^0 Z^0 (\Pi(X^0) - \Pi^{\text{sol}}). \quad (23) $$

In the above expressions, $\kappa^1$ and $\kappa^0$ are nonnegative functions of the densities and the electrical potential at the interfaces, while $Z^0$ and $Z^1$ are nondecreasing functions of the jump of grand potential at the interfaces satisfying $Z^0(0) = Z^1(0) = 0$. The parameters $\Pi^{\text{sol}}$ and $\Pi^{\text{met}}$ stand for the grand potential in the solution and in the metal respectively and are assumed to be given.

The fact that the displacement of the interfaces is driven by the grand potential seems to be new in the context of corrosion modeling. However, such a fact was suggested in the seminal book (Pimpinelli and Villain, 1999), and used in relatively close contexts by (Cermelli and Jabbour, 2005; Li et al., 2009), but in the absence of electrical effects.

## 2.2 Onsager's reciprocal relation and free energy dissipation

The model we propose has been designed in order to ensure the decay along time of some thermodynamical potential, namely the Gibbs free energy of the full system made of the solution, the oxide and the metal.

Define $G_{ox}(t)$ the Gibbs free energy of the oxide at time $t$ as the sum of a chemical contribution $A_{ox}(t)$ and an electrical contribution $\mathcal{E}_{ox}(t)$ respectively defined by

$$ \mathcal{A}_{ox} = \int_{X^0}^{X^1} \left( A_{hl} + \sum_{i=1}^{3} A_i(u_i) \right) dx \quad (24) $$

and

$$ \mathcal{E}_{ox} = \frac{\lambda^2}{2} \left( \int_{X^0}^{X^1} |\partial_x \Psi|^2 dx + \frac{1}{\gamma^0} (|\Psi^0|^2 - |\Delta \Psi_{pzc}^0|^2) + \frac{1}{\gamma^1} (|\Psi^1|^2 - |V - \Delta \Psi_{pzc}^1|^2) \right) \quad (25) $$

with $\Psi$ deduced from the $u_i$ thanks to the Poisson equation (1)–(3) detailed previously. Then one shows that $\frac{\delta A_{ox}}{\delta u_i} = \mu_i$ and $\frac{\delta \mathcal{E}_{ox}}{\delta u_i} = z_i \Psi$, so that $\frac{\delta G_{ox}}{\delta u_i} = \xi_i$ owing to (16). Here $\frac{\delta G_{ox}}{\delta u_i}$ stands for the first variation of $G_{ox}$ with respect to $u_i$. Therefore, the time evolution of the Gibbs free energy writes

$$ \frac{d}{dt} G_{ox} = \sum_{i=1}^{3} \int_{X^0}^{X^1} \xi_i \partial_t u_i dx. \quad (26) $$

After a few calculations to be detailed in a forthcoming article building on the model equations described in the previous section, one gets that

<figure><img src="image_12.png" /></figure>
$$ \frac{d}{dt} G_{ox} = \int_{X^0}^{X^1} \sum_{i=1}^{3} G_i \partial_x \xi_i dx - \sum_{i=1}^{3} (F_i^1 \xi_i^1 - F_i^0 \xi_i^0) - R_{PB} \dot{X}^1 \pi^1 + \pi^0 v_d^0. \quad (27) $$

Each contribution in the integral of the right-hand side is nonnegative because of (17) and (18). The other terms correspond to exchange terms with the solution and with the metal. Their sign is not clear, motivating the introduction of an augmented free energy accounting for the solution and the metal. We set

$$ G_{tot}(t) = G_{ox}(t) + \int_{0}^{t} \sum_{i=1}^{3} [F_{i}^{1}(\tau)\xi_{i}^{\text{met}} - F_{i}^{0}(\tau)\xi_{i}^{\text{sol}}]d\tau + (X^{1}(t) - X^{0}(0))[R_{PB}\Pi^{\text{met}} + (1 - R_{PB})\pi^{\text{sol}}] - (X^{0}(t) - X^{0}(0))\Pi^{\text{sol}}. $$

Then it follows from the calculations detailed above that

$$ \frac{d}{dt} G_{tot} = \int_{X^0}^{X^1} \sum_{i=1}^{3} G_i \partial_x \xi_i dx + \sum_{p \in \{0,1\}} \sum_{i \in \{1,3\}} F_i^p \Delta \xi_i^p + F_2^1 \Delta \xi_2^1 + F_2^{0,PR} \Delta \xi_2^{0,PR} + F_2^{0,FR} \Delta \xi_2^{0,FR} + R_{PB} \dot{X}^1 \Delta \Pi^1 + v_d^0 \Delta \Pi^0. $$

In the above right-hand side, for $i = 1$ and $i = 3$, we have set $\Delta\xi_i^0 = \xi_i(X^0) - \xi_i^{\text{sol}}$ with $\xi_i^{\text{sol}} = \log(\frac{m_i^0}{k_i^0})$ being the electrochemical potential of species $i = 1$ and $i = 3$ in the solution. Similarly, for $i \in \{1,2,3\}$, we have set $\Delta\xi_i^1 = \xi_i^{\text{met}} - \xi_i(X^1)$, with $\xi_i^{\text{met}} = \log(\frac{k_i^1}{m_i^1}) + z_i V$. For the fluxes of electrons across the solution/oxide interface, we have to introduce two electrochemical potentials $\xi_2^{0,PR}$ and $\xi_2^{0,FR}$, but the reasoning is similar. Then as highlighted for instance in (Cancès et al., 2023a), the structure of the Butler-Volmer laws (4)-(5) and (7)-(10) is such that all the contributions corresponding to fluxes across the interfaces in the above right-hand side are nonnegative.

To establish the decay of the augmented free energy along time, it only remains to check that the two last terms related to the motion of the interfaces yield nonnegative contributions. This stems from the fact that the functions $Z^p$ are equal to 0 at 0 and increasing.

To conclude this section, let us remark that the decay of the Gibbs free energy for our model is a consequence of the particular mathematical structure of our model. Indeed, our model has been derived as the generalized gradient flow of the augmented Gibbs free energy (Mielke, 2011; Peletier, 2014; Peletier et al., 2022).

# 3. A reduced model with fixed interfaces

## 3.1 Model simplification

In the spirit of previous contributions of the literature – see for instance (Chainais-Hillairet and Lacroix-Violet, 2014, 2012) – we have focused on the simplified model where oxygen vacancies are neglected. Thanks to this simplifications, interfaces can be assumed to be fixed, and after rescaling the width of the oxide layer can be assumed to be constant along time and equal to 1. This simplification may sound strong, but the reduced model allows to give good insights for the production rate of dihydrogen $H_2$. Moreover, it can be used for the calibration of the full model and its numerous parameters.

## 3.2 Mathematical analysis

In the reduced model we proposed in (Cancès et al., 2023a), the unknowns are $u_1$ and $u_2$, as well as the self-consistent electrostatic potential $\Psi$ solving the Poisson equation. There is no further need of the oxygen vacancy density $u_3$, neither on the grand potential $\Pi$, but the other ingredients remain the same: Butler-Volmer law for the fluxes across the interface, nonlinear convection-diffusion equation for the

<figure><img src="image_13.png" /></figure>
ferric cations and linear convection-diffusion equation for the electrons. Furthermore, the system can still be interpreted as the generalized gradient flow of some augmented Gibbs free energy. This ensures in particular some stability for the model on which we build to establish the global in time existence of a solution to the system in (Cancès et al., 2023a).

It is worth noticing that we do not impose restrictions on the parameters, in opposition to what was necessary in the previous work (Chainais-Hillairet and Lacroix-Violet, 2014). Moreover, we also establish the physically relevant bounds $0 < u_1 < \bar{u}_1$, $0 < u_2$ and $0 < u_3 < \bar{u}_3$ among other technical estimates we do not detail here, cf. (Cancès et al., 2023a) for interested readers. These bounds are all consequences of the control of the free energy and of its dissipation rate, following the methodology of (Gajewski and Gröger, 1989; Jüngel, 2015).

## 3.3 Numerical approximation

As a step towards the development of a finite volume scheme for the full model, we addressed the reduced model. The preservation of the consistency of our model with the second principle of thermodynamics at the discrete level has been a priority when designing the numerical method. This led us to consider two-point flux approximation (TPFA) finite volumes in space and a backward in time Euler time discretization. For the fluxes in the drift-diffusion equations for the conservation of the charge carriers, we either use so-called Scharfetter-Gummel (SG) fluxes for the linear equation corresponding to electrons, or a new extension of the squareroot approximation (SQRA) fluxes to the case of nonlinear mobilities for ferric cations (Cancès and Venel, 2023). We also developed and analysed mathematically in (Cancès et al., 2023b) a new variant of the SG scheme which shares the main features with the new SQRA scheme, that are the preservation of the physical bounds and the decay of some discrete free energy functional.

Even though this would not lead to particular problem, we did not prove rigorously so far the convergence of our scheme towards the global in time solution to the reduced model exhibited in (Cancès et al., 2023a). We rather focused on the simulation of the reduced model thanks to a prototype code developed in Python, and to its comparison with the model obtained from the original DPCM (Bataillon et al., 2010) after similar simplifications consisting in neglecting the oxygen vacancies. We refer to (Cancès et al., 2023a) for the values assigned to the parameters in our tests.

<figure><img src="image_14.png" /><figcaption>Figure 3 – Comparison of the current-voltage characteristics corresponding to the new thermodynamically consistent reduced model (vDPCM, orange solide line) and to the original reduced model (DPCM, dashed blue line) for pH = 7 (left), pH = 8,5 (center) and pH = 10 (right).</figcaption></figure>

It appears that from a macroscopic point of view, our new reduced model (referred to as vDPCM in what follows) gives similar results as the reduced model building on the original model (still referred to as DPCM in the figures below). This can for instance be observed on the i-v curves depicted on Figure 3. However, it clearly appears on Figure 4 that the original model does not preserve the physical bound $u_i < \bar{u}_i$ on the density of ferric cations, in opposition to our new model.
<figure><img src="image_15.png" /></figure>

Figure 4 – Profiles of the scaled densities of cations (top), electrons (middle) and the electrostatic potential (bottom, in Volts, the physical unit) at the steady state for two values of the applied potential V, that are -0.4 Volts (left) and 0.3 Volts (right).

# 4. Conclusion and prospects

## 4.1 Conclusion

We proposed a new model which can be seen as an update of the model proposed in (Bataillon et al., 2012, 2010) to accurately represent the magnetite layer protecting the iron canister to be stored in geological repositories. Our new model is designed to be compatible with the second principle of thermodynamics. More precisely, we built the model as a generalized gradient flow so that Onsager reciprocal relation is automatically satisfied.

We performed in (Cancès et al., 2023a) the mathematical and numerical analysis of a simplified model where oxygen vacancies have been neglected, so that the geometry of the domain does not evolve along time. The new model we propose provides close results from a macroscopic point of view on the few cases we assessed, but significant differences are observed in terms of the repartition of the charge carriers in the oxide layer. In particular, our model preserves the physical range for the ferric cation occupation rate, in opposition to the simplified model (Chainais-Hillairet and Lacroix-Violet, 2014) built from the original model (Bataillon et al., 2010).
New contributions to the theory of finite volumes have been needed to properly handle the nonlinear mobilities in the convection-diffusion of cations. We proposed in (Cancès et al., 2023b; Cancès and Venel, 2023) some schemes that are free energy diminishing while keeping a small computational cost.

Let us finally mention the work (Merlet et al., 2022) where another simplification is addressed. A reduced model accounting only for oxygen vacancies (but with moving interfaces) is proposed and analysed thanks to tools from calculus of variations and Wasserstein gradient flow theory {Citation}.

## 4.2 Prospects

We still have to develop an efficient numerical tool to simulate the full model sketched in section 2.1. On top of the non-trivial treatment of the time discretization, the discretization of the grand potential $\Pi$ requires a particular attention. The highly coupled and nonlinear nature of the problem makes the effective resolution difficult as we enforce strong constraints on the computational cost and the robustness of the tool so that it can be coupled to THMC platforms without significantly increasing the computational time.

Among the other challenges we will have to face in the near future, the calibration of our model has to be done carefully. Indeed, the calibration procedure that was employed to fix the parameters in the original DPCM model cannot be used directly and must be adapted. In particular, the choice of the functions $Z^0, Z^1$ and $\kappa^0, \kappa^1$ governing the motion of the interfaces is left open so far. Advanced discussions with electrochemists and specialists of material science will be required here.

Once our Python prototype code will be finalized, then we plan to use it to train a surrogate model to get a fast solver to be used in the context of uncertainty quantification and sensitivity analysis.

## Code source

We make use of a private Python prototype code which is not mature enough to be published on public repositories. The code will be made available when finalized and calibrated.
## References

Bataillon, C., Bouchon, F., Chainais-Hillairet, C., Desgranges, C., Hoarau, E., Martin, F., Perrin, S., Tupin, M., Talandier, J., 2010. Corrosion modelling of iron based alloy in nuclear waste repository. Electrochimica Acta 55, 4451–4467. https://doi.org/10.1016/j.electacta.2010.02.087

Bataillon, C., Bouchon, F., Chainais-Hillairet, C., Fuhrmann, J., Hoarau, E., Touzani, R., 2012. Numerical methods for the simulation of a corrosion model with moving oxide layer. Journal of Computational Physics 231, 6213–6231. https://doi.org/10.1016/j.jcp.2012.06.005

Breden, M., Chainais-Hillairet, C., Zurek, A., 2021. Existence of traveling wave solutions for the Diffusion Poisson Coupled Model: a computer-assisted proof. ESAIM: M2AN 55, 1669–1697. https://doi.org/10.1051/m2an/2021037

Cancès, C., Chainais-Hillairet, C., Merlet, B., Raimondi, F., Venel, J., 2023a. Mathematical analysis of a thermodynamically consistent reduced model for iron corrosion. Z. Angew. Math. Phys. 74, 96. https://doi.org/10.1007/s00033-023-01970-6

Cancès, C., Herda, M., Massimini, A., 2023b. Finite volumes for a generalized Poisson-Nernst-Planck system with cross-diffusion and size exclusion.

Cancès, C., Venel, J., 2023. On the square-root approximation finite volume scheme for nonlinear drift-diffusion equations. Comptes Rendus. Mathématique 361, 535–558. https://doi.org/10.5802/crmath.421

Cermelli, P., Jabbour, M., 2005. Multispecies epitaxial growth on vicinal surfaces with chemical reactions and diffusion. Proc. R. Soc. A. 461, 3483–3504. https://doi.org/10.1098/rspa.2005.1495

Chainais-Hillairet, C., Gallouët, T.O., 2016. Study of a pseudo-stationary state for a corrosion model: Existence and numerical approximation. Nonlinear Analysis: Real World Applications 31, 38–56. https://doi.org/10.1016/j.nonrwa.2016.01.010

Chainais-Hillairet, C., Lacroix-Violet, I., 2014. On the existence of solutions for a drift-diffusion system arising in corrosion modeling. DCDS-B 20, 77–92. https://doi.org/10.3934/dcdsb.2015.20.77

Chainais-Hillairet, C., Lacroix-Violet, I., 2012. The existence of solutions to a corrosion model. Applied Mathematics Letters 25, 1784–1789. https://doi.org/10.1016/j.aml.2012.02.012

Gajewski, H., Gröger, K., 1989. Semiconductor Equations for variable Mobilities Based on Boltzmann Statistics or Fermi-Dirac Statistics. Mathematische Nachrichten 140, 7–36. https://doi.org/10.1002/mana.19891400102

Jüngel, A., 2015. The boundedness-by-entropy method for cross-diffusion systems. Nonlinearity 28, 1963–2001. https://doi.org/10.1088/0951-7715/28/6/1963

Li, B., Lowengrub, J., Rätz, A., Voigt, A., 2009. Geometric evolution laws for thin crystalline films: modeling and numerics. Commun. Comput. Phys 6, 433–482.

Merlet, B., Venel, J., Zurek, A., 2022. Analysis of a one dimensional energy dissipating free boundary model with nonlinear boundary conditions. Existence of global weak solutions. HAL:03888607. https://hal.science/hal-03888607

Mielke, A., 2011. A gradient structure for reaction-diffusion systems and for energy-drift-diffusion systems. Nonlinearity 24, 1329–1346. https://doi.org/10.1088/0951-7715/24/4/016

Peletier, M.A., 2014. Variational Modelling: Energies, gradient flows, and large deviations. arXiv:1402.1990 [math-ph].

Peletier, M.A., Rossi, R., Savaré, G., Tse, O., 2022. Jump processes as generalized gradient flows. Calc. Var. 61, 33. https://doi.org/10.1007/s00526-021-02130-2

Pimpinelli, A., Villain, J., 1999. Physics of Crystal Growth, Physics of Crystal Growth.