Melpomeni Dimopoulou
email: dimopoulou@i3s.unice.fr

Eva Gil

San Antonio

Marc Antonini

A quaternary code mapping resistant to the sequencing noise for DNA image coding

The exponential growth in the generation of digital information creates a big challenge for data storage given the capacity limitations of conventional storage devices. Recent works have proposed DNA as a means of digital data storage proposing a novel solution for long-term storage. Although having many advantages, DNA storage is a challenging topic due to the error-prone process of DNA sequencing (reading). To deal with this error most existing works focus on the introduction of error-correction methods. However, most of those methods introduce important redundancy without promising full error correction for the widely used sequencing method using the Nanopore sequencer. This work focuses on noise resistance rather than error-correction proposing a new algorithm for optimally assigning VQ indices to DNA codewords while reducing the visual impact of substitution errors that are caused during sequencing.

I. INTRODUCTION

Data explosion is leading to an increasing demand of storage resources which is inconsistent with the capacity and longevity limits of traditional storage units. Interestingly enough, studies estimate that 90% of the internet data has been only generated in the last two years while 80% of this information consists the so called "cold" data [START_REF] Shacklett | Unstructured data: A cheat sheet[END_REF] which is very rarely or even never accessed but still needs to be safely stored for security and compliance reasons. The archiving of cold data requires tons of hardware which should be replaced every 5 to 10 years to ensure reliability. Every year millions of data centers worldwide are wasting huge amounts of energy and money for the migration of data into new back-up drives in order to protect all this information. Therefore, it is clear that the need of finding new ways for storing data is crucial and to this end recent studies have proposed the storage of data in form of DNA strands as a potential alternative to the current methods. The biological properties of the DNA molecule allow an efficient data storage for hundreds of years without loss of information and without additional costs to maintain its correctness. An interesting example of DNA's incredible longevity is a woolly mammoth which had been trapped into permafrost whose DNA has been successfully decoded after 40.000 years. DNA data storage is clearly a very promising multidisciplinary solution which consists of both data science algorithms and some delicate biological processes. More This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 863320.

precisely, the DNA data storage can be described by the following main steps. To begin with, the digital data has to be transformed into a quaternary encoding of the symbols A, T, C and G which denote the four different building blocks (nucleotides) of DNA. This sequence is then synthesized into DNA by the biological process of DNA synthesis. This process is performed in vitro and is error free as long as the DNA strands to be synthesized are shorter than 300 nucleotides (nts). Consequently, before synthesis the encoded quaternary strand needs to be cut into small chunks (oligos) and formatted to contain special headers that denote the position of each oligo in the initial long sequence. The synthesized oligos are then stored into some safe containers (capsules) which prevent contacts with water and oxygen and promise long term preservation. To retrieve the stored information in the decoding part of the workflow, the stored DNA strands will be read in vitro using the biological process of DNA sequencing. The sequencers are special machines which can read the DNA content of the capsules providing to the output the quaternary data. Finally, the quaternary sequence is transformed back to its initial representation to retrieve the stored information. The biggest challenge of this workflow is imposed by the process of DNA sequencing as the sequencers can introduce errors in the decoding process. This error can be reduced if the encoding performed in the first step of the workflow respects the following rules:

• No homopolymers: A nucleotide should not be repeated more than three times in a row. • G,C content: % G,C ≤ % A,T • No pattern repetition: Repetition of short patterns should be avoided.

Consequently, the encoding algorithm should be carefully selected to provide an encoding that respects all the above rules. In this paper, we present an extension of our previously proposed encoding [START_REF] Dimopoulou | Efficient Storage of Images onto DNA Using Vector Quantization[END_REF] by applying a more sophisticated mapping function which is resistant to substitution errors which can be created during the sequencing. This special mapping allows the reduction of the total distortion in the decoded image in the case that it is corrupted by errors. In section II we present some interesting works proposed by the state of the art. In section III we describe the algorithm which is used for constructing a robust quaternary code while in section IV we explain the mapping algorithm that is used for the encoding into a quaternary sequence. More precisely in IV-A we introduce some basic notions to facilitate understanding of the following section IV-B where we present the proposed algorithm for optimally mapping input quantization vectors to the quaternary codewords. In section VI we present the methods used in our experiments and we demonstrate our experimental results. Finally in section VII we discuss about the conclusions and future works.

II. STATE OF THE ART

As discussed in the previous section DNA data storage is a very promising yet very challenging process as the decoding of DNA is error-prone and thus the encoding is constrained by important restrictions and should provide a code which is robust to errors. In the state of the art there have been proposed some very interesting algorithms for encoding an input stream into a robust quaternary representation ( [START_REF] Goldman | Towards practical, high-capacity, lowmaintenance information storage in synthesized DNA[END_REF], [START_REF] Blawat | Forward error correction for DNA data storage[END_REF], [START_REF] Grass | Robust chemical preservation of digital information on DNA in silica with errorcorrecting codes[END_REF], [START_REF] Erlich | DNA fountain enables a robust and efficient storage architecture[END_REF], [START_REF] Yazdi | A rewritable, random-access DNA-based storage system[END_REF], [START_REF] Dimopoulou | A biologically constrained encoding solution for long-term storage of images onto synthetic DNA[END_REF]). In our latest work in [START_REF] Dimopoulou | Efficient Storage of Images onto DNA Using Vector Quantization[END_REF] we proposed the use of Vector Quantization to further improve the efficiency of the compression. In this work we propose an extension of our solution in [START_REF] Dimopoulou | Efficient Storage of Images onto DNA Using Vector Quantization[END_REF] by adding a special mapping algorithm which is resistant to noise. This mapping is strongly inspired by the works in [START_REF] De Marca | An algorithm for assigning binary indices to the codevectors of a multi-dimensional quantizer[END_REF] and [START_REF] Farvardin | A study of vector quantization for noisy channels[END_REF] where the authors proposed an algorithm for assigning binary words to vectors of a multidimensional quantizer in such a way so that in case of an error the decoded vector would be close to the correct one improving this way the quality of the decoded signal. In this paper we extend this algorithm to a constrained quaternary dictionary to be applied for DNA coding.

III. CREATING THE QUATERNARY CODE

For the creation of the quaternary code we will use the algorithm proposed in our previous work in [START_REF] Dimopoulou | A biologically constrained encoding solution for long-term storage of images onto synthetic DNA[END_REF]. The main idea of this algorithm lies on the construction of words using some pair-symbols, the concatenation of which can provide viable words. As viable words we consider the words which when assembled in an encoded strand they provide a quaternary sequence which respects all the constraints imposed by the process of DNA sequencing and will therefore be robust to sequencing noise. The codewords of the code are build using permutations of the elements from the two following dictionaries:

• D 1 = {AT, AC, AG, T A, T C, T G, CA, CT, GA, GT } • D 2 = {A, T, C, G}
More precisely, words of an even length l are constructed using all possible permutations of l 2 pair-symbols from dictionary D 1 . Equivalently, words of an odd length l are constructed using all possible permutations using l-1 2 pair-symbols from D 1 to represent the first l -1 nucleotides and a single symbol from D 2 to represent the last nucleotide of the words. This method for building quaternary words ensures that the codewords produced contain an appropriate percentage of G and C which does not exceed 40-60% while also avoiding homopolymers. It is clear that due to the biological restrictions imposed, some quaternary words will be omitted from the produced DNA code as they are considered as non-viable for being more prone to sequencing errors. Consequently, the constructed code will be less efficient than a complete quaternary code but it is more robust to distortion. Another important asset of this algorithm is the low complexity while it provides a valid fixed length code that can be applied to any type of input data in contrast to most of the existing state of the art algorithms that can only be applied to binary.

IV. A CONTROLLED CODE-MAPPING RESISTANT TO SEQUENCING NOISE

A. Introduction to the proposed mapping

Sequencing of DNA strands is an error-prone process and can create insertions, deletions and substitutions of nucleotides causing important noise in the decoded data. Illumina and Nanopore are the two most widely used types of sequencers. While Illumina is the most reliable solution providing a low error-rate in the decoding, scientists are turning towards the more recent Nanopore machine due to its low cost, fast throughput and user-friendly, small size. Unfortunately this type of sequencer introduces more noise in the decoded sequence compared to the Illumina, creating the need for finding new ways to deal with errors in the decoded sequence. Existing studies have proposed some interesting methods for error correction by adding redundancy in the encoded sequence, while always using robust encoding to reduce the possibility of sequencing noise. However, the high error-rate introduced by sequencers can't easily be completely eliminated. Thus, in this study we take a very first step in proposing an algorithm which provides a solution resistant to sequencing noise. In other words, this algorithm aims to reduce the impact of the remaining sequencing errors. In this first attempt, we will only focus on the noise of substitutions.

In [START_REF] De Marca | An algorithm for assigning binary indices to the codevectors of a multi-dimensional quantizer[END_REF] DeMarca et al. have proposed an algorithm for assigning binary words to codevectors of a multi-dimensional quantizer in such a way so to be resistant to channel noise. Inspired by this idea and under the assumption that the noise of substitutions introduced by sequencers can be modeled as the one introduced by noisy channels of transmission we extend the algorithm proposed in [START_REF] De Marca | An algorithm for assigning binary indices to the codevectors of a multi-dimensional quantizer[END_REF] which is applied to binary codewords to an application which is using a quaternary code. Before describing the proposed mapping algorithm it is important to introduce some basic notions and definitions. Let V = {v 1 , v 2 , . . . , v K } be a set of -dimensional vectors v k . To encode those K vectors, a code W = {w 1 , w 2 , . . . , w K } of K quaternary codewords is necessary. Since this study focuses on DNA data storage this code W is constructed using the dictionaries D 1 and D 2 providing a constrained encoding as explained in section III. Consequently, the size K of the code W is restricted to specific values as the words should be created according to the constraints imposed by the process of DNA coding. Consequently K ∈ {k 1 , k 2 , . . . } where:

k 1 = 10 and k i+1 = 4k i , i: odd 10k i-1 , i: even
Given the two sets V and W, one can define a mapping function M : V → W for assigning codewords from W to codevectors of V. Using this function one can encode a sequence of input blocks (codevectors) to a quaternary representation of A, T, C and G.

Assuming a substitution error in the encoded quaternary stream, some correct codeword w c will be transformed to an erroneous codeword w e and under the hypothesis that the error rate produced by the sequencer will be reasonably small, the Hamming distance between the correct and the erroneous codeword will be d H (w c , w e ) = 1. Therefore, for each correct codeword w c of length l there are different erroneous codewords of distance 1. The set of all possible erroneous codewords can be, according to coding theory, represented in the Hamming Space as a sphere H(w c ) of radius 1 the center of which is the correct codeword w c . An example of such a sphere is depicted in figure 1. In the case where the code W contained all the possible arrangements of A, T, C, G, there would be 4 l different spheres the cardinality of which would be |H(w c )| = 3l. However, as explained in section III, since the code used for DNA coding excludes some words which can't be viable due to the encoding constraints, in this work we consider K different spheres with varying cardinality. In other words, some codewords w e that would normally belong to some sphere of center w c might be omitted due to the fact that they do not respect the rules of DNA coding. As a result, a substitution can cause two different possible types of error:

• Decodable error: The substitution transforms a correct codeword w c to an erroneous word w e which exists in the constrained code C proposed in section III and therefore w e ∈ H(w c ). Decoding will then provide an erroneous vector v e instead of the correctly decoded vector v c . In the case that the Euclidean distance d E (v c , v e ) is small, the produced error will not significantly affect the visual quality of the decoded image. To this end we propose a special mapping algorithm which allows the assignment of codeword indices to vectors such that the possible errors will lead to a minimum distortion. This algorithm is analytically described in section IV-B. • Undecodable error: The substitution transforms a correct codeword w c to an erroneous codeword w e which does not exist in the constrained code C. In this case decoding is not possible and thus the application of some error correction is necessary to allow decoding. The applied correction techniques are further described in section V.

B. The proposed controlled mapping algorithm

The main goal of this algorithm is to map codewords w q from a code W which differ from each other at exactly one position (Hamming distance of 1) to codevectors v q from a vector set V which are close in terms of Euclidean distance. The objective of this mapping lies in the fact that in case of an error during sequencing and assuming that the sequencing noise is small enough, a correct codeword will be transformed to another one which will have a small Hamming distance with the correct one. Thus, if those two codewords are assigned to input vectors which are close, the error will not significantly affect the image quality. To further define which vectors can be considered as close vectors we introduce a set S(v) which contains the closest vectors to the vector v in terms of Euclidean distance. The central idea of the proposed mapping is to assign codewords of the same sphere to vectors which belong to the same neighborhood as shown in figure 3. However, such an assignment is not possible for every neighborhood S and thus it is necessary to perform this assignment according to some priority. To this end, we define an empirical function F (v) introduced by [START_REF] De Marca | An algorithm for assigning binary indices to the codevectors of a multi-dimensional quantizer[END_REF] for a vector v as:

F (v) = p(v) α β (v)
where p(v) is the probability of v in the input sequence, and

α(v) = j|vj ∈S(v) d E (v j , v)
with β ≥ 0 a trade-off parameter. Therefore vectors with a higher value for F are considered to belong in a denser neighborhood and should consequently get higher priority to be assigned to the same sphere of words. The algorithm can be very roughly described by the following parts:

• For each codeword w q : Create a sphere containing the B q codewords which have a Hamming distance of 1 compared to w q . Define B = max i (B i ), i ∈ {1, . . . , L}.

• For each input vector v q :

-Find a set S(v q ) of B neighboring vectors v n which are the closest to v q in terms of Euclidean distance

d E (v q , v n ). -Compute the following empirical function F (v q ).
• Use algorithm 1 to progressively perform assignment of vectors v q to codewords w q such that vectors with a bigger F (v q ) as well as their neighboring vectors v n ∈ S(v q ) will be assigned to the same sphere of codewords whenever possible as depicted in figure 1. If this is not possible assignment is performed such that vectors are assigned to codewords which have a small Hamming distance from the codewords already assigned to their neighboring vectors. • Optimization of the first assignment:

-Exchange the previously mapped codewords between each pair of vectors.

-For each exchange check if the average distortion has decreased. If true keep this change, else keep the initial state of mapping.

V. PROPOSED DECODING OF UNDECODABLE WORDS

As discussed in section IV-A in the case in which a substitution error creates an undecodable word it is necessary to employ some error correction to allow decoding. In this section we will describe two possible methods for correcting undecodable words. The first method which we will call Simple Correction Decoding (SCD) is working according to the two following assumptions.

• A substitution error will produce an erroneous word which will not differ from the original one in more than one nucleotide. • Using the mapping algorithm proposed in the previous section it is probable that codewords that are close in terms of Hamming distance are assigned to codevectors which are close in terms of Euclidean distance. Let us assume a word w c which is corrupted by a substitution error producing an undecodable erroneous codeword w e . Given the first assumption, d H (w c , w e ) = 1. To allow decoding, w e should be first corrected to a decodable word w d . The idea behind SCD is to correct w e to some word w d we chose as most appropriate the codeword w d ∈ W which is assigned to the median vector-index of the input image. This choice has two possible outcomes. Either w d will be equal to w c creating no visual distortion in the decoded image, or

w d = w c and d H (w c , w d ) = {1, 2}.
By making use of one extra assumption which indicates that since the input data is an image there can be correlations between neighboring elements we propose a second more sophisticated decoding, which we will call Advanced Correction Decoding (ACD). ACD works in two decoding cycles. In the first cycle it performs a first decoding by omitting the undecodable words. Then in a second round the algorithm decodes the remaining undecodable words using the following steps: -Define as V H the set of vectors v Hj to which have been assigned the words in H(w e ).

-Define W n as the set of neighboring words w ni around w e with 0 ≤ i ≤ 8 as depicted in figure 3 -Define set V n as the set of vectors v ni with 0 ≤ |V n | ≤ 8 to which the codewords in W n have been assigned.

-Define set S = W n ∩ H(w e ) = {w 1 , w 2 , . . . , w z } • If |S| = 1, w e ← w 1 ∈ S • If |S| > 1 -Define f (w z ) as the frequency of a word w z , w z ∈ S -w e ← w z ∈ S such that w z = arg max z f (w z ) • if |S| = ∅ -Compute D(v Hj ) = |Vn | i=1 d E (v H j ,vn i )f (M -1 (vn i )) |V H j | , for v Hi ∈ V H , ∀v Hj ∈ V H -w e ← w d where w d = M(v d ) with v d := arg min v H j (D(v Hj ))

VI. EXPERIMENTAL RESULTS

In our study we performed multiple comparisons in order to prove the efficiency of the proposed mapping algorithm. For each comparison we ran 10 different realisations of random error added on the same input image and have plot the improvement in PSNR in function of the error rate. The points in each plot correspond to the mean value of the 10 different realisations of noise. As explained in IV-B, the main assumption behind the proposed mapping algorithm is the restriction of the errors to one error per quaternary word. Thus, for the noise addition, we used a uniform random distribution to select a percentage of words w c from the encoded sequence and we introduced an error of one nucleotide in each w c . To begin with, to set an upper bound in the performance of the controlled mapping we tested the improvement of PSNR compared to the non-controlled case by adding only decodable errors. More specifically, in this experiment we added random substitutions making sure that the created erroneous words exist in the constrained codebook. Figure 4a, shows the improvement of PSNR (∆PSNR) for different parameters of K and for VQ in function of the introduced error rate. As observed in this figure, the controlled mapping can improve the PSNR by at least 3 and at most 7.5 dB in the optimal case where no corrections are needed. The addition of only decodable errors might not be a realistic case in practice but this plot reveals the true potential of the proposed mapping algorithm. This is because in the case of an undecodable error the correct codeword will be transformed to an erroneous one w e with d H (w c , w e ) = 1. This second more realistic experiment is presented in figure 4b which depicts the improvement of PSNR in function of different error rates comparing the controlled and non-controlled mapping using SCD for correcting the undecodable words. As expected the undecodable errors decreases the improvement of PSNR providing a best case of ∆PSNR= 3dB. To improve those results and since SCD is a simple decoding where the undecodable word is corrected to the codeword mapped to the median vector-index of the quantized input image, we also tested the case of ACD. The result is shown in figure 4c which depicts the improvement in PSNR comparing the controlled mapping using ACD to the non-controlled mapping using SCD. This case reveals the total improvement which can be achieved providing a ∆PSNR betwenn 2.5 and 5dB. On the right column of figure 4 we also present the visual improvement for each of the previous cases on the input image which was used in our experiments. The visual results are presented for an error rate of 3%, which is the percentage of substitutions introduced by the sequencer of the Nanopore (according to a study in [START_REF] Zeng | Causalcall: Nanopore basecalling using a temporal convolutional network[END_REF]), and for the values of K and which provided the best PSNR improvement in each case.

VII. CONCLUSION

In this work we have introduced a new noise resistant controlled mapping for optimally assigning DNA quaternary words to vectors obtained by VQ. This new mapping has been tested on a constrained codebook according to the needs of DNA data storage. The obtained results exhibit a very promising increase of PSNR providing a noticeable visual improvement. An interesting future step of this work would be extending this algorithm to treat other kinds of errors such as insertions and deletions.

Set V = {v1, v2, . . . , vK } of codevectors v k , |V| = K Set W = {w1,

1 Fig. 1 :

 11 Fig.1: Example of a Hamming sphere. The cross-elements denote non viable words that would belong in the Hamming sphere but are omitted due to the constrained quaternary code.

Fig. 2 : 1 Fig. 3 :

 213 Fig. 2: Mapping vectors to codewords

F (v)} 4 :Fig. 4 :

 44 Fig.4: Left column: PSNR improvement in function of the error rate for different cases of VQ parameters K and with β = 0. This value of β was reported in[START_REF] De Marca | An algorithm for assigning binary indices to the codevectors of a multi-dimensional quantizer[END_REF] to perform optimally for a non-uniform input probability distribution. The three different plots correspond to different experiments as explained in the legends a, b and c. Right column: Visual results for each of the cases a, b and c for the values of K and that provided the best performance. The selected error rate for those images is 3% accordingly to the percentage of substitutions caused by the Nanopore sequencer[START_REF] Zeng | Causalcall: Nanopore basecalling using a temporal convolutional network[END_REF].

978-1-7281-9320-5/20/$31.00 c 2020 European Union


Simone Guiso

Dung Di Caprio

Jacques De Lamare
email: jacques.de-lamare@cea.fr

Benoît Gwinner

Intergranular corrosion: comparison between experiments and cellular automata

Keywords: intergranular corrosion, cellular automata, modelling, material ageing, stainless steel

In this paper, a 3D cellular automata approach of intergranular corrosion in presented. A model is proposed, along with experimental results both qualitatively (groove morphology, groove angle distribution) and quantitatively (mass loss, corrosion rate, corrosion front evolution). Simulations have been performed using a new asynchronous algorithm: the choice of the most adequate grid is also discussed. The model investigates the statistical variation of the corrosion front and provides useful information to make predictions on material ageing.

Introduction

In France, the PUREX (Plutonium and Uranium Rening by EXtraction) process is used to reprocess spent fuel in a nuclear fuel cycle [START_REF] Fauvet | Corrosion issues in nuclear fuel reprocessing plants[END_REF]. The process involves the treatment of the spent fuel in hot and concentrated nitric acid for liquid-liquid extraction of the recoverable material (uranium and plutonium).

In some particular conditions, the interaction between the nitric media and the structural stainless steel (SS) may lead to intergranular corrosion (IGC) phenomena [START_REF] Fauvet | Corrosion issues in nuclear fuel reprocessing plants[END_REF][START_REF] Andreoletti | Industrial feedback of intercycle uranyl nitrate concentration in stainless steel 304L evaporator[END_REF][START_REF] Ueno | Corrosion phenomenon of stainless steel in boiling nitric acid solution using largescale mock-up of reduced pressurized evaporator[END_REF].

Understanding IGC is therefore crucial. This particular type of corrosion promotes grain detachment and leads to corrosion kinetics larger than expected.

IGC study is a challenge that involves both experimental campaigns and simulation modelling.

Experimentally, IGC corrosion kinetics is mainly studied using weight-loss (WL) techniques and electrochemical measurements. However, experiments are intrinsically limited in duration and cannot always provide sucient information to predict material behaviour over long time periods. Modelling can then provide an important support.

The way IGC nds its path along grain boundaries is essential to understand IGC corrosion kinetics and may lead to dierent characteristic morphologies in relation to grain detachment. Thus, on the modelling side, dierent paths have been recently investigated, where focus is on the evolution of the corrosion morphology. In brick-wall models [START_REF] Zhao | Generalized model for 38 IGC growth in aluminum alloys[END_REF], the IGC penetration is described as a function of the probability of choosing a dened path in a "wall structure". Monte-Carlo simulations were performed by Lim et al. [START_REF] Lim | Overview of intergranular corrosion mechanisms, phenomenological observations, and modeling of AA5083[END_REF] to estimate the average IGC penetration depth at precise time steps in comparison with experimental tests. A stochastic reduced order model (SROM) was instead proposed by Sarkar et al. [START_REF] Sarkar | Stochastic reduced order models for uncertainty quantication of intergranular corrosion rates[END_REF] to quantify the uncertainty on IGC rates. Results demonstrate that this approach is more suitable than MC method for this purpose. Despite all these achievements, the IGC kinetics is yet to be fully understood.

In particular, more investigations should be addressed to the morphology characterization of the IGC grooves and to their progression inside the material, which is dicult to access experimentally. A complete description of the phenomenon would be an important tool for the industry in strong need of predictive models for material ageing. The objective of this study is to start lling this gap by presenting a predictive tool for material ageing that includes a ne description of the material morphology evolution. In this paper, the study is focused on "non-sensitized" SS. A 3D stochastic CA model is presented and compared with some experimental tests already published [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF][START_REF] Tcharkhtchi-Gillard | Kinetics of the oxidation of stainless steel in hot and concentrated nitric acid in the passive and transpassive domains[END_REF].

The paper is organized as follows: a brief background to the IGC study 2 is presented in Section 2, followed by a description of the experimental tests in Section 3, and of the CA model in Section 4, while results are discussed in Section 5. The conclusive remarks are then presented in Section 6. Tables of symbols, abbreviations and references are provided at the end.

Background

From the experimental point of view, IGC has been studied in terms of WL measures. Specimens are usually exposed to a corrosive solution and extracted at precise time steps [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF][START_REF] Whillock | Intergranular corrosion testing of austenitic stainless steels in nitric acid solutions[END_REF]. The goal is to obtain information on mass loss and corrosion rate as a function of time.

In parallel, electrochemical tests are also performed to evaluate the risk of IGC.

More generally, electrochemical measurements are classically used to characterize dierent aspects of the SS corrosion. For example, the determination of the pitting potential denes the sensitivity to this type of corrosion [START_REF]Standard Test Method for Conducting Cyclic Potentiodynamic Polarization Measurements for Localized Corrosion Susceptibility of Iron-, Nickel-, or Cobalt-Based Alloys[END_REF], double loop electrochemical potentiokinetic reaction tests inquire about the degree of sensitization of the SS [START_REF]Standard Test Method for Electrochemical Reactivation (EPR) for Detecting Sensitization of AISI Type 304 and 304L Stainless Steels[END_REF], electrochemical impedance spectroscopy experiments give information on the limiting corrosion processes [START_REF] Orazem | Electrochemical Impedance Spectroscopy[END_REF]. Among those, linear voltammetry is often performed in the case of IGC to investigate whether the SS remains passive (where the SS dissolves slowly and uniformly) or it is polarized in its transpassive state (where IGC occurs) in its chemical environment (see Figure 1) [START_REF] Ningshen | Corrosion degradation of AISI type 304l stainless steel for application in nuclear reprocessing plant[END_REF][START_REF] Diaz | Corrosion performance of AISI-309 exposed to molten salts V 2 O 5 -Na 2 SO 4 at 700 • C applying EIS and R p electrochemical techniques[END_REF].

In order to resist corrosion, the SSs used for structural material are screened in view of their corrosion resistance, optimized and qualied before their use. In particular, austenitic SSs with low carbon, high chromium and controlled silicon and phosphorus contents are preferred [START_REF] Armijo | Intergranular corrosion of nonsensitized austenitic stainless steels[END_REF]. A stable passive layer protects all these materials so that corrosion, in normal conditions, is slow and uniform. However, too oxidizing conditions (e.g. nitric media containing oxidizing ions at high temperature like Vanadium (V)) may promote the deterioration of this passive layer and generate IGC phenomena. In this case, the cathodic reaction (reduction reaction) is fast and can lead to a shift of the corrosion potential to Normally, grain boundaries (GB) represent barriers that confer strength to metallic materials. However, in some environmental conditions, they are the main source of weakness and failures. Grain boundary degradation includes corrosion, cracking, embrittlement and fractures. According to literature [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF][START_REF] Callister | Materials science and engineering: an introduction[END_REF] two cases of IGC can be distinguished.

The rst one includes the "sensitized" SSs, which present a relatively high carbon content. When SSs are subjected to long heat treatments at high temperatures (for instance, between 500 and 800 • C), Cr and C species may diuse through the material and precipitate as chromium carbides (Cr 23 C 6 ) at grain boundaries. As chromium slowly diuses from the grain matrix, a poor-inchromium-content region is generated adjacent to Cr 23 C 6 precipitates. In this region, the chromium content is not sucient to provide a stable passive lm to protect the SS. As a consequence, the SS is exposed to IGC attack [START_REF] Callister | Materials science and engineering: an introduction[END_REF]. To prevent these phenomena, dierent countermeasures have been considered. In the industrial eld, for example, titanium (Ti) or niobium (Nb) are added to the material. Their presence in the matrix promotes the formation of Ti x C or Nb x C precipitates and prevents the formation of chromium carbides. Thus, the material maintains its passive layer -rich in chromium -and IGC attacks are avoided [START_REF] Singh | Research on the diusion of ti-based stabilizer through ux coated GTAW to reduce carbide precipitation in SS304 weldment[END_REF][START_REF] Pinnow | Corrosion resistance of stainless steel weldments[END_REF].

As an alternative, the carbon content in SS may be reduced in order for it to be lower than 0.03 wt.% [START_REF] Callister | Materials science and engineering: an introduction[END_REF]. In this case, the carbon is soluble in the austenitic structure and it remains uniformly distributed in the crystal structure. The SS is then called "non-sensitized" and chromium carbides formation is avoided during high temperature treatments. However, IGC can still occur in some specic conditions. In the context of investigations on the corrosion behavior of a SS in a spent fuel reprocessing plant [START_REF] Fauvet | Corrosion issues in nuclear fuel reprocessing plants[END_REF], experimental tests showed that the material can still be aected by IGC in presence of highly oxidizing media. The degradation of a non-sensitized SS is characterized by the formation of triangular grooves at grain boundaries (as represented in Figure 2). The progression Figure 2: Cross section of a "non-sensitized" 310L SS suering IGC [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF]. of the grooves in the material bulk may cause grain detachment with a certain periodicity, as shown in Figure 3.

On the modelling side, as already stated in the introduction, dierent paths Figure 3: Scanning Electron Microscope observation of detached grains of a "non-sensitized" 304L SS [START_REF] Tcharkhtchi-Gillard | Kinetics of the oxidation of stainless steel in hot and concentrated nitric acid in the passive and transpassive domains[END_REF].

have been recently investigated, among those a promising method is to model IGC phenomena through the cellular automata (CA) approach [START_REF] Di Caprio | 3D cellular automata simulations of intra and intergranular corrosion[END_REF]. It consists of a dynamical system where both space and time are discretized [START_REF] Pérez-Brokate | Cyclic voltammetry simulations with cellular automata[END_REF][START_REF] Chopard | Cellular Automata Modeling of Physical Systems[END_REF]. Each site (named cell) is identied by a state (dened from a nite list of possible states) and its spatial coordinates. The dynamics of the system is then regulated by a list of transition rules, which take into account the interaction between neighbouring cells (generally rst-neighbour cells) [START_REF] Cordoba-Torres | Non-diusive spatial segregation of surface reactants in corrosion simulations[END_REF].

Recently, dierent studies have been performed within this approach. Lishchuk et al. [START_REF] Lishchuk | A cellular automaton model for predicting intergranular corrosion[END_REF] improved their previous model [START_REF] Lishchuk | Model-based prognosis for intergranular corrosion[END_REF] by developing a 3D stochastic CA model. A similar result was obtained by Jahns et al. [START_REF] Jahns | Modeling of intergranular oxidation by the cellular automata approach[END_REF] who, on the other hand, developed a 2D CA model to study the intergranular oxidation phenomenon.

Igarashi et al. [START_REF] Igarashi | Simulations of intergranular corrosion feature for stainless steel using cellular automata method[END_REF][START_REF] Ueno | Study of the eect of phosphorus concentration on intergranular corrosion of stainless steel in boiling nitric acid solution[END_REF] also focused their studies on the development of a 2D CA model to study IGC in stainless steels, while Di Caprio et al. [START_REF] Di Caprio | 3D cellular automata simulations of intra and intergranular corrosion[END_REF] developed a 3D stochastic CA model to study IGC, by improving the previous 2D model [START_REF] Taleb | Numerical simulation of the eect of grain size on corrosion processes: Surface roughness oscillation and cluster detachment[END_REF].

Gwinner et al. [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF][START_REF] Gwinner | A stereological approach for measuring the groove angles of intergranular corrosion[END_REF] proposed a phenomenological geometrical model to estimate the long-term corrosion rate of "non-sensitized" SS suering IGC. The model assumptions are based on experimental observations: grooves grow symmetrically along the grain boundaries, groove angles are constant with time surface. Similar observations were obtained by Beaunier et al. [START_REF] Beaunier | A kinetical model for the electrochemical grooving of grain boundaries[END_REF], which proposed to describe IGC with two principal parameters:

V gb which represents the grain boundary dissolution rate V s which represents the grain dissolution rate

The dierent approaches are summarized in the following Table 1 and2, along with their main characteristics.

Determination of the experimental corrosion velocities

Figure 4: Schematic representation of the corrosion rates (V gb and V s , where

V gb > V s ), the groove depth (H * ), the groove angle (α) and the groove inclina- tion (β) [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF].

Experimental tests were performed by Gwinner et al. [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF] on an AISI 310L stainless steel. It presents a low carbon content and was rstly thermally treated (homogenization treatment and quenching) to avoid precipitate formation that would have brought a type of IGC proper to sensitized SS (see Section 2). The SS was exposed to a solution of nitric acid (HNO 3 ), containing an oxidizing ion (VO + 2 ) at boiling temperature, that was renewed after each corrosion period of 160 h. SS specimens were periodically removed and examined. The experimental equivalent thickness loss (ETL exp , in µm) was evaluated with the following relation

ETL exp (t) = m(t) ρS 10 4 (1)
where m(t) is the cumulated mass loss (in g) at time t, ρ is the SS density (8 g.cm -3 ) and S is the initial surface of the polished sample (in cm 2 ) that is in contact with the corrosive agent. The instantaneous experimental corrosion rate CR exp is then dened as

CR exp t + ∆t 2 = ETL exp (t + ∆t) -ETL exp (t) ∆t ( 2
)
where ∆t is the duration of a corrosion period between two mass loss measurements.

Besides the experimental tests, a semi-empirical model was also developed to estimate the long term IGC evolution for the SS [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF]. According to this model, a triangular groove is represented as in Figure 4.

The model is based on experimental observations of the mass loss kinetics.

According to experimental data, the corrosion rate (V corr ) linearly evolves from an initial value to a steady-state condition. The initial value corresponds to the metal dissolution at the initial state (t=0), whose value is V s . The steadystate value corresponds to the corrosion rate at which the IGC kinetics is fully controlled by grain dropping: in this case, the corrosion rate is equal to V * gb . As discussed in [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF], the transition time from the linear trend to the steady-state condition is estimated (in good agreement with the experimental observations) to be at t = 1.5

D V * gb
, where D is the average grain dimension.

Following these assumptions, the corrosion rate (V corr ) can be expressed as:

t ≤ 1.5 D V * gb , V corr = 2 3 V * gb D V * gb -V s t + V s (3) t ≥ 1.5 D V * gb , V corr = V * gb (4)
V s and V * gb are then obtained by tting the rst 4 points of the experimental corrosion rate kinetics with the linear trend of the semi-empirical model.

V gb is then determined by the following formula:

V gb = V gb * cos(β) ( 5
)
where β is dened in Figure 4 and averaged over dierent grooves [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF].

In [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF], V s = 92µm/y, V * gb = 1287µm/y and β = 35 • have been obtained.

3D Cellular automata model

The 3D stochastic CA model will hereafter be presented. The methodology belongs to the so-called "grey-box" class of approaches as it stands between "black-box" and "white-box" models: the former are phenomenological approaches where the model response is directly related to input data but no interpretation nor mechanistic description is introduced; the latter are based on a detailed understanding of physical laws and processes, but are too complex to be implemented. The CA model presented in this paper adopts a simple stochastic mesoscopic level of description. The main features of the model are those who characterize the CA approach: a lattice grid, a list of possible cell states, the corrosion rules. The model is then related to "real" space-time quantities by providing space and time scales. Finally, the construction of the Voronoï diagram, used to model the granular microstructure of the material, is detailed.

Cell states, corrosion rules and probabilities

The CA model is made of a grid of cells that can acquire three possible chemical states: IGN, which represents the grain boundaries, GRN that represents the grains and SOL that denes the corrosive agent.

At each iteration i, CA cells can change their state according to the transition rules. In the case of the corrosion process, IGN and GRN cells change their states to SOL according to two eective probabilities, P ign for grain boundaries and P grn for grains, as presented in Figure 5. When a SOL cell comes into contact with a GRN cell (or IGN cell), a random number δ is drawn and compared to a probability P grn (or P ign ) associated to the cell according to the following equations

P grn = P grn N b N b,max (6a)
P ign = P ign N b N b,max (6b)
where P grn and P ign are the intrinsic corrosion probabilities for the GRN and

IGN cells. The ratio N b N b,max
is a correction related to the cell environment, If δ < P grn (similarly δ < P ign ), the cell will change its state to SOL, i.e. the cell is "corroded", otherwise it keeps its state (and the cell is not "corroded"). In this CA model at mesoscopic scale, the goal of Eq. (6a) and (6b) is to represent the fact that the more SOL cells surround a GRN (or an IGN) cell, the faster corrosion is. The lifetime of a cell is assumed to be directly proportional to the number of surrounding SOL cells, considering the fact that dissolution is proportional to the surface in contact with the solution. This feature is new in such a type of model.

In order to compare the model to experimental data, the choice of both intrinsic corrosion probabilities P grn and P ign is crucial. A preliminary study of pure grain dissolution is rstly performed to assess the functional relation between probabilities and corrosion velocity in the simulation. A system with a single 3D block of GRN cells surrounded by three layers of SOL cells is adopted.

The precise number of SOL layers is in fact immaterial as long as all GRN top plane cells are in contact with at least a SOL cell. In this study, the P grn probability is then varied. Following the results in [START_REF] Guiso | Inuence of the grid cell geometry on 3D cellular automata behavior in intergranular corrosion[END_REF], a hexagonal close-packed (HCP) grid is chosen for this preliminary study.

Figure 6: Grain dissolution velocity v stat grn and chemical roughness λ stat in stationary regime as function of the grain corrosion probability P grn .

In Figure 6, the results of the study are presented in terms of corrosion velocity v stat grn (P grn ) (expressed as the number of GRN cells dissolved per iterations per cell in the horizontal plane) and chemical roughness λ stat (P grn ), as a function of P grn in stationary conditions. The chemical roughness, also known as rugosity in some domains, is the number of GRN cells in contact with at least one SOL cell, normalized with respect to the horizontal surface of the system. This quantity is of interest because it is directly related to the reactivity of the surface, as it gives information on the kinetics. This contrasts with the standard physical roughness (generally dened as the root mean square of the deviation of the prole with respect to the mean height) that is more a morphological information. Note also that the physical roughness increases as a power law with time; a stationary regime is reached only due to nite size eects which limit the power law. Corrosion velocity and chemical roughness are averaged between the iteration n, which indicates the beginning of the stationary regime, and iteration N , in such away that N -n is a suciently large number of iterations for the average. Their expressions read as

v stat grn (P grn ) = 1 N XY N grn (n) -N grn (N ) N -n (7)
λ

(P grn ) = 1 N XY N i=n N grn-sol (i) N -n (8)
where N grn (i) is the number of GRN cells that are still present in the system after i iterations and N grn-sol (i) is the number of GRN cells that are in contact with at least one SOL cell at the iteration i, N XY is the total number of cells that are present in a horizontal XY -plane.

Figure 6 shows that the velocity is a linear function of P grn and that the chemical roughness remains constant with P grn .

Figure 7 shows another characterisation of the surface: the percentage histogram of the number of SOL neighbours for the GRN cells on the surface.

Such information is relevant as the probabilities (Eq. ( 6a) and (6b)) depend on the number of neighbours to a cell. The statistical distribution of all SOL neighbours is shown to remain identical whatever the value of P grn . The results above are important in order to set the relations between simulated and experimental quantities and allow the comparison between the two.

The linearity of the corrosion velocity with P grn indicates that changing P grn is equivalent to scaling the real time duration of a simulation iteration. This means that, in order to represent a given corrosion rate, one can rst associate a simulation iteration i to a given time duration and then tune the value of P grn to obtain the experimental corrosion rate. However, owing to linearity, there is not one sole choice. The same process can be represented by increasing the duration of an iteration (for instance multiplying it by a given coecient) and correlatively decreasing P grn (that is dividing it by the same coecient). Note that, alternatively, it is possible to choose P grn and the iteration time step in accordance to the corrosion rate. In all cases, the system is the same and the quantity of dissolved material per real time units remains unchanged. If the described system is identical, the chemical roughness and the SOL neighbours distribution remain unchanged, independently of P grn : which is indeed what has been proved in the simulations.

Note that the linearity between corrosion rate and probabilities is not evident, as not all algorithms provide such property. In this paper, this property is a consequence of the asynchronous algorithm (Figure 5), where cell states are updated in the order as they come along during simulations. The standard synchronous algorithm (where all lattice cells updating is performed simultaneously), that was chosen in [START_REF] Guiso | Inuence of the grid cell geometry on 3D cellular automata behavior in intergranular corrosion[END_REF], does not exhibit linearity between corrosion velocity and P grn for all values. A synchronous algorithm was also used in a previous IGC study in [START_REF] Di Caprio | 3D cellular automata simulations of intra and intergranular corrosion[END_REF].

An important consequence of the linearity of the corrosion rate with P grn , is that it is possible to model a given v grn v ign velocity ratio using an identical P grn P ign probability ratio.

In conclusion, the two CA corrosion probabilities (P grn and P ign ) are chosen as follows. To achieve faster simulations, the highest probability is associated to the fastest process. Adopting a smaller value is less ecient as some iterations do not lead to a IGN dissolution event. The IGN corrosion probability P ign is set equal to 1. Then, P grn is set in order to have P grn P ign ratio equal to the experimental ratio v grn v ign , which in this case corresponds to V s V gb ≈ 0.0581 (see 3).

Grid

The choice of an appropriate grid is important to best describe the behaviour of a material. The specic kinetics of IGN corrosion is related to the crystallographic misorientation of adjacent grains and to the chemical composition of these regions. GB are schematically represented as planes and do not have a preferential spatial orientation. Therefore any realistic model must statistically reproduce GB corrosion velocities, which are isotropic and independent on the planar GB spatial orientation. Dealing with CA approach means that simulations are performed on discrete lattices, which therefore are intrinsically anisotropic. A recent study [START_REF] Guiso | Inuence of the grid cell geometry on 3D cellular automata behavior in intergranular corrosion[END_REF] aimed at nding the best suited lattice to get as close as possible to an isotropic behaviour. It showed that the HCP grid is the closest to an "isotropic" grid in view of IGC corrosion along grain boundaries. However, this study was performed using a synchronous algorithm, which has been here substituted by an asynchronous algorithm. The purpose of this section is to investigate whether an identical statement can be valuable with the present asynchronous algorithm. data generated with four dierent random seeds. For a grain boundary thickness of two cells, the average deviation ε is equal to -2.05%, which indicates that simulation velocity is slightly lower than the theoretical one. The standard deviation is around 1.58%. Figure 10 shows that the deviation δ stands within 2%, except for particular angles at θ = 90 • and ϕ = 0 [90 • ], which are characterized by a higher velocity ("red points" in the Figure). To understand these specicities, the layout of the grain boundaries for these particular angles was investigated in detail. It was found that the velocity increase is caused by a thickness eect, which occurs when the IGN plane is parallel to a grid face.

In these geometries, despite an integer grain boundary thickness of two, the grain boundary of the HCP lattice appears to be constituted of a slab with an alternation of two and three cells thickness. In the model, the IGN dissolution probability (Eq. 6b) depends on the number of solution neighbours. A thickness larger than two in these directions causes an increase of the velocity. This eect was investigated by varying the GB slab thickness.

In Figure 11, two dierent orientations have been investigated, (θ, ϕ) = (90 • , 0 • ) and (90 • , 30 • ). These orientations of the grain boundary correspond to the most extreme cases representing respectively the highest and the lowest velocities with respect to the average. For both orientations, the velocity is an increasing function of thickness. Such an outcome is expected since the probability (Eq. 6b) is an increasing function of the number of neighbours, which is larger when thickness increases. However, the velocity dierence between the two orientations decreases when the thickness increases. This can be understood by the relative dierence in thickness between ϕ = 0 • and ϕ = 30 • that decreases when the slab thickness increases. Note that results are averaged here over 50 dierent simulations: as a consequence, the error bars are small.

In conclusion, the relative small variations between the simulated and the theoretical velocity seen in Figure 10 for asynchronous algorithm conrm the result that was obtained for the synchronous algorithm in [START_REF] Guiso | Inuence of the grid cell geometry on 3D cellular automata behavior in intergranular corrosion[END_REF]: the HCP grid is then chosen for the comparison with the experimental data.

Space and time, CA to experiments conversion factors

A summary of the main parameters for the experimental setup and the corresponding CA simulations is presented in Table 3. The grain size within the CA model is set by the number of initial seeds N seed , independently on the grain dimension in the "real" experimental case: in this case N seed = 104. This value corresponds to a balance between having, within the simulation box, both suciently large grains and a sucient number of grains in order to describe the corrosion progression in the material. The characteristic dimension of a grain is determined by the number of seeds with the following equation:

D sim = 3 N xyz N seeds = 272 (10)
where N xyz is the total number of cells that are present in the system.

To compare experimental data and simulation results, it is necessary to associate the real space and time scales with CA dimensionless quantities. In the CA approach, these scales correspond respectively to the cell distance and to the duration of one iteration. At this point, two conversion factors need to be introduced.

The rst coecient is a space scaling factor A (in µm/cell), which can be calculated by the ratio between the experimental average grain dimension (in µm) [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF] and the CA average grain dimension (in number of cells). In this case,

A = 84 272 ≈ 0.309 µm/cell.
Once the spatial scale is dened, the second coecient is the time scaling factor B (in hour/iteration) that enables to convert the duration of a simulation the corrosion experiments). The coecient A is also present because the ratio between experimental velocity and CA velocity involves a spatial ratio, with the space expressed once in µm and once in cells. Note that the specic value for B equal to 1 is not meaningful; it is a coincidence arisen from the experimental conditions and the chosen units.

Voronoï tessellation

The material microstructure has been reproduced through a Voronoï diagram [START_REF] Okabe | Locational optimization problems solved through voronoi diagrams[END_REF][START_REF] Voronoi | Nouvelles applications des paramètres continus la théorie des formes quadratiques. deuxième mémoire, recherches sur les parallélloèdres primitifs[END_REF]. It consists of a partitioning of the 3D grid into regions originated by a set of starting points, called seeds. To build the Voronoï diagram, a spherical region is progressively grown around each seed along all directions: when two regions meet each other, a border is determined along the mediator plane.

Then, each seed keeps growing in all directions until all grid cells become part of a region including a seed or a border. Note that periodic boundary conditions are applied by replicating the entire seed spatial distribution to complete the Voronoï structure at boundaries in all directions. Once the Voronoï tessellation is complete, states are assigned to cells: those which belong to borders assume the IGN state, the others the GRN state. This construction leads to grain boundaries with a precise two cells thickness, reducing the above-mentioned thickness velocity issue.

The initial grid is then obtained by assigning the SOL state to the three uppermost XY -layers: during dissolution, the solution path is conducted from the top of the grid to the bottom. Note that the top layer grain size distribution is dierent than the one in the bulk because the grains on the surface have been

cut. An example of the initial grid is represented in Figure 12.

Figure 12: Representation of the 3D initial grid of IGC model.

Computational setup

The simulation code has been written in C language, in NVIDIA's CUDA environment, and simulations have been performed using NVIDIA Tesla K80 cards in Dell PowerEdge C4130 servers with Intel Xeon E5-2640 processors. A single simulation run takes from 24 to 48 hours, depending on the values that are given to P grn and P ign (the higher the values, the lower the simulation time).

To obtain a reasonable statistics, results are averaged over 50 dierent Voronoï diagrams. A comparison between a section of the experimental reference test specimen and a section of the 3D CA model is presented in Figure 13. Both structures of grains are similar, although some dierences can be observed. Some grain boundaries are not clearly visible on the experimental image due to a lower reactivity of these grain boundaries during the metallographic etching process.

Results

Preliminary observations on the morphology

Some crystal twinning appears on the experimental image, which are not modelled by CA. A majority of triple grain boundaries forming 120 • angles suggests the structure is stabilized. Figure 14 shows 3D perspectives and 2D sections of CA simulations during corrosion at t = 400, 800, 1200, 1600, 2000 h. The groove progression inside the material is visible as well as the consequent grain detachment. Good qualitative agreement with experimental observations is found (cf. Figure 2). Some detached grains are presented in Figure 15, in comparison with Figure 3.

Figure 15: 3D view of some detached grains at t=1600 h.

Figure 16 shows the distribution of the α angles, the α-projection observed along some 2D cross sections of the system for both experimental tests and CA simulations. Two "real" times, with the correspondent iterations, have been considered: the rst one, t=641h, in the groove progression regime, while the second one, t=1602h, includes the groove progression and the grain detachment phenomena. CA α angle distribution is an average over all α angles that are observed in 2 equally spaced cross sections for 10 dierent simulations. In each section, 7 to 8 grooves were generally found.

Results conrm that, in both approaches, the distribution does not have a signicant evolution in time. In the simulations, a narrow peak was found between 5 • and 15 • , similar to the experimental observations [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF]. A comparison of the main results is presented in Table 4: the dierence between experimental and CA average values of α is within the standard deviation. A distribution tail for larger α is observed. It is due to the fact that α angles are obtained from 2D sections and are therefore larger than the corresponding α 3D angles, as discussed in [START_REF] Gwinner | A stereological approach for measuring the groove angles of intergranular corrosion[END_REF]. Figure 17 shows the comparison between simulations and experiments, in terms of ETL and CR. CA simulation results have been scaled with the abovementioned space and time factors. The equivalent thickness loss has been di-vided in two components: the ETL due to the material dissolution (represented in black) and the one due to grain dropping (dierence between the red squared points and the black curve). This distinction is not possible in the experiments.

In order to match the denitions that were used in [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF], the following corresponding denitions are considered

ETL sim (t) = N mat (0) -N mat (t) N xy A (11)
where N mat (t) is the sum of the number of IGN and GRN cells at time t, N XY is the number of grid cells in a XY -plane and A is the space scaling factor.

Note that, as the horizontal axis are now plotted in real time, iteration i can be replaced with time t, as B = 1.

The corrosion rate is dened by

CR sim t + ∆t 2 = ETL sim (t + ∆t) -ETL sim (t) ∆t × 8760 B ( 12
)
where ∆t is chosen equal to 160 (a time interval large enough to average corrosion rates), the factor 8760 is introduced again to have standard corrosion rates in µm/y. The 3D CA model shows a good agreement with the experimental tests. At the beginning, the corrosion rate grows linearly with time while the equivalent thickness loss squarely evolves with it. Then, a transition regime is observed between 600 h and 1000 h. After t=1000 h, the corrosion rate becomes constant and the weight loss increases linearly with time. This behaviour is related to the morphology of the system. Firstly, grooves are formed and progress into the material, causing a progressive grain detachment. Then, grains start dropping with a certain periodicity (depending on the size). On average, when the rst grain layer is detached, a stationary corrosion rate is reached.

Evolution of the surface morphology

The considerations on the mass loss discussed in the previous paragraph give only a mean information of the degradation of the SS in terms of mean thickness loss. However, taking into account the localized character of IGC, the corrosion is deeper in some locations of the SS like at the bottom of the grooves, as illustrated in Figure 2 and Figure 14. As these points dene the residual available thickness of SS, it is fundamental to have information on their time-dependent position. However, information on the deepest corrosion is experimentally only accessible through analysis of 2D cross sections (Figure 2)

which destroy the material. The evolution of the relative corrosion surface gives a rst set of information on this aspect. It is dened as the ratio of the corroded surface at time t to the corresponding initial planar surface. In practice it can be calculated from 2D cross section pictures, as the ratio of the contours of the corroded surface to their horizontal projection. A comparison between numerical and experimental results is given in Fig. 18. The numerical results were obtained by averaging over four equally spaced cross sections and ten Voronoï structures.

Both curves present a similar trend. The relative surface evolution starts with a transitory regime up to around t=500 h; then the stationary regime starts. Once again, the trend of this parameter reects the two corrosion regimes: groove progression (transitory) and grain detachment (stationary).

CA results are comparable with experimental ones. The dierences between the two curves is within the error bars. At a given corrosion time, the spatial distribution of the deepest points of the grooves corresponds to the spatial distribution of the grain boundary cells that are in contact with the solution (IGN-SOL). An example of an IGN-SOL interface is presented in Figure 19. It should be noted that the synchronous algorithm adopted in [START_REF] Di Caprio | 3D cellular automata simulations of intra and intergranular corrosion[END_REF] for P ign = 1 showed a homogeneous dissolution of IGN with the lowest points at identical height. The present asynchronous algorithm is more satisfactory in this case, since it shows a distribution of such points similar to experimental observations. CA simulations consider, at each (x,y) coordinate, the lowest IGN cell that is in contact with a SOL cell, while experiments consider only the lowest points of the grain boundary grooves. At low depth, a narrow peak is present; then, as the depth increases, the number of non-independent events increases and the distribution spreads. When the rst grain characteristic dimension is reached on average in terms of depth and grain detachment starts, a small hump appears aside the main distribution peak.

The evolution of the IGN-SOL interface depth as a function of time is presented in Figure 21. The three curves represent respectively the minimum (in red), the average (in blue) and the maximum (in black) depths, estimated every 25 simulations (which correspond to 25 hours). The average and the maximum interface depths present a linear trend. The minimum interface depth presents instead a linear trend up to 500 h, then it switches to another linear trend with a slightly higher slope. This particular behaviour is again linked to the two regimes that characterize the IGC phenomenon. These data have a high industrial interest, since they give important information on when a material may reach a critical failure point in terms of corrosion. The time evolution of the standard deviation (SD) of the IGN-SOL interface depth distribution is shown in Figure 22. At rst, it seems to follow a linear trend; then it switches to a square-root-of-time trend. The two trends correspond to the two morphological regimes that characterize the IGC, as seen in Figure 20. The linear regime matches with the groove progression inside the material, while the second one reects the combination of the previous eect together with grain detachment. Due to the limited dimensions of the simulated systems, the grain detachment events have a large impact on the corrosion progression and explain the oscillations at t > 1000 h.

CA modelling has been used as a so called grey box approach. Like for all simulations, it is necessary to introduce in the simulation the real space and time scales. For the space, the relation between grain size and number of cells is chosen. There is no unique choice: a "best choice" is made to suciently balance a detailed description (large lattice) and the simulation time (not too large lattice). The time scale is then derived from a velocity ratio using the stationary IGC velocity as reliable physical quantity, as it is obtained as an average over several points. Finally, the grain corrosion velocity is taken from the experimental study. It should be noted that, as a grey box approach, the method cannot predict corrosion rates, but it uses coecients provided experimentally. However, by only using three input numbers (two of them (A and B) are necessary to give a meaning to the simulation) several results (in the form of trends or graphs) are obtained. The main results are listed below.

The model reproduces several IGC features that are compared to experimental results of [START_REF] Gwinner | Towards a reliable determination of the intergranular corrosion rate of austenitic stainless steel in oxidizing media[END_REF]: the Voronoï tessellation reasonably accounts for the material morphology;

the approach shows the characteristic grooves expected for "nonsensitized" SS IGC corrosion; the simulation provides the angular distribution of the grooves; the ETL and corrosion velocity curves are well reproduced and exhibit the two characteristic corrosion kinetic regimes; the relative surface morphology is also reproduced.

New features are obtained: information about grain dropping and pure granular dissolution have been separated; lowest corrosion point depth distribution histogram is traced; supplementary information of minimum and maximum of these points and their standard deviation have been obtained.

Note that some features are dicult to obtain experimentally as they concern the bulk material and require the physical destruction of the sample as they cannot be performed in situ.

Conclusion

In this paper, a 3D cellular automata model is proposed to investigate the IGC phenomenon. The model provides realistic morphological characteristics of IGC like groove shaped corrosion, grain detachment. The model can be considered as a grey box approach. The approach feeds on experimental parameters and combines them with mechanistic features within the CA approach. Here, corrosion velocities are used to parametrize the corrosion probabilities of the model. In contrast with [START_REF] Di Caprio | 3D cellular automata simulations of intra and intergranular corrosion[END_REF][START_REF] Guiso | Inuence of the grid cell geometry on 3D cellular automata behavior in intergranular corrosion[END_REF], a new asynchronous algorithm has been used which ensures linearity between CA parameters and real corrosion kinetics. This is useful to treat IGC modelling which involves dierent corrosion kinetics for granular and intergranular material.

Concerning the quantitative comparison, a particular attention has been spent to validate the choice of the CA grid. The model can thus provide quantitative kinetics of the ETL and angular distribution of the grooves. The model also gives new insights on IGC as it is able to separate the contributions from granular and intergranular corrosion and gives the statistical spatial distribution of the penetration depths of the corrosion. The precise knowledge of this distribution is of critical interest for managing intergranular corrosion of SS in industry. These features help characterizing the material dissolution and assess system reliability predictions.

This article shows how to pass from experimental to simulation quantities.

The model is simple and draws input parameters directly from phenomenological corrosion velocities. However, ner mechanistic descriptions can be considered in the future. For instance, the corrosion probabilities can be based on intrinsic quantities like energetic barriers and allow taking into account temperature eects. A detailed description of the environment can also be examined. In this model, a simplied eect of the strength of the solution is introduced by connecting the corrosion probability to the amount of solution surrounding a corroded cell. In the future, it will also be possible to account for diusion and concentration eects of corrosive compounds in solution.

Figure 1 :

 1 Figure 1: Simplied scheme of the dierent electrochemical domains.

Figure 5 :

 5 Figure 5: Schematic representation of the corrosion probabilistic transition rule.

Figure 7 :

 7 Figure 7: Distribution of the number N b of SOL neighbours to GRN cells during grain dissolution for dierent values of P grn .

Figure 8 :

 8 Figure 8: Schematic 3D view of a planar GB (in blue) dened by its polar and azimuthal angles (θ, ϕ) and its normal n [31].

Figure 9 :

 9 Figure 9: Schematic representation of the HCP grid, with the representation of the vertical cell stacking.

Figure 10 :

 10 Figure 10: Simulated corrosion velocity results for HCP grid, expressed in terms of δ(θ, ϕ) (in percentage colored scale).

Figure 11 :

 11 Figure 11: Simulated corrosion velocity as a function of the grain boundary thickness for two congurations: (θ, ϕ) = (90 • , 0 • ) in blue and (θ, ϕ) = (90 • , 30 • ) in red. Results are averaged over 50 dierent simulations.

Figure 13 :

 13 Figure 13: Microstructure before corrosion of AISI 310L stainless steel [7] (left) in comparison with a cut view of the Voronoï microstructure (right).

Figure 14 :

 14 Figure 14: 3D views (left) and 2D sections for the planes indicated in red (on top) of the system at times t = 400, 800, 1200, 1600, 2000 h. Detached grains are removed for visibility.

Figure 16 :

 16 Figure 16: Comparison of experimental (in pink and red) and CA simulated (light-blue and blue) histogram distributions of groove angles at dierent time steps, t=641 h and t=1602 h, and at corresponding CA iterations. For each set of bars, α values are referred to the α interval dened by the values in the abscissa inside which they are represented.

Figure 17 :

 17 Figure 17: Weight loss comparison between CA simulations (in red) and experimental tests (in blue). The black line represents the portion of equivalent thickness loss that comes from the material dissolution.

Figure 18 :

 18 Figure 18: Comparison of the relative surface with experimental results in blue and simulated results in red. The error bars are estimated at 2σ/ √ Z, where σ is the standard deviation and Z is the number of tests.

5. 4 .

 4 Evolution of the distribution of the deepest penetrations

Figure 19 :

 19 Figure 19: Representation of the grain boundary -solution interface during corrosion.

Figure 20 :

 20 Figure 20: IGN-SOL spatial distribution as a function of the depth for experimental results (red shades) et CA simulations (blue shades). Lines are provided as guides.

Figure 21 :

 21 Figure 21: Representation of the average IGN-SOL interface depth as a function of time. Results are calculated every 25 simulation steps.

Figure 22 :

 22 Figure 22: Representation of the standard deviation of the IGN-SOL spatial distribution as a function of time.

Table 1 :

 1 Summary of the dierent IGC models with their main features

	Model	Approach	Results
			Minimum IGC path length as a function
			of time
	Brick-wall [4]	2D, probabilistic
			Comparison with experimental results
			Comparison with experimental results
	Monte-Carlo [5]	3D, probabilistic	IGC damage depth as a function of the
			degree of sensitiveness and of time
			Comparison SROM/Monte-Carlo simula-
			tions
	SROM [6]	1D, mathematical and probabilistic model	Eect of the uncertainty in an-odic/cathodic sites on the current
			density
			Eect of the domain length on the current
			Groove angular description
	Geometrical model [7, 29]	2D/3D, deterministic	Comparison with experimental results
			Corrosion depth and corrosion rate

Table 2 :

 2 Summary of the dierent CA IGC models with their main features

	Reference	Approach	CA grid / cell geometry	Results
				Comparison with experimental results
	[23]	2D, probabilistic	Cubic / cubic 2D: 1216 cells/grain	Corrosion depth, path length
			3D: 448000 cells/grain
				Comparison with experimental results
	[24]	3D, probabilistic	Cubic / cubic 3D: 448000 cells/grain	Corrosion depth, path length
				Corrosion path
	[25]	2D, probabilistic	Square / Square	Comparison with experimental results
			Hexagonal grains	(only microstructure)
				Comparison with experimental results
	[26, 27]	2D, probabilistic	Square / square Grain size randomly	Corrosion width and depth, corrosion rate
			distributed
				Corrosion depth and corrosion rate
	[19]	3D, probabilistic	Cubic / cubic Grain size randomly	Roughness study
			distributed
				Evolution of the corrosion front as a
				function of dierent spatial arrange-
			Square / square	ments (square, brick-wall, polygonal)
	[28]	2D, probabilistic	Various types of grains	Surface roughness prole
			and spatial distributions
				Corrosion rate

Table 3 :

 3 Comparison between experimental and simulation main parameters.

	Parameter	Experimental tests [7]	CA simulations
	Dimension	30 × 20 × 1.5 mm 3	1280 × 1280 × 1280 cells
	Average grain lateral dimension	84 µm	272 cells
	Stationary vertical corrosion velocity	1206 ± 152 µm/y	0.446 ± 0.003 cell/iteration
	Corrosion velocities/	v gb = 1582.3 ± 187 µm/y,	P ign = 1,
	probabilities	v s = 92 ± 46 µm/y	P grn ≈ 0.0581
	step into "real" hours. It is calculated as the ratio between the stationary
	CA corrosion velocity and the experimental one. In this case, B = 8760 ×
	0.446 × 0.309 1206	≈ 1.

0 hour/iteration, where an extra factor 8760 accounts for the change of time unit from years (used in corrosion velocity) to hours (used in

Table 4 :

 4 Comparison between experiments and simulations in terms of average alpha and standard deviation at two dierent time steps

	Parameter	Experimental tests [7] t = 641 h t = 1602 h i = 641 i = 1602 CA simulations
	Average alpha	11.6 •	11.5 •	10.9 •	9.9 •
	Standard deviation	6.2 •	5.1 •	3.8 •	3.6 •
	5.2. Mass loss kinetics

Table 5 :

 5 List of symbols S initial surface of polished sample in cm 2 D exp , D sim characteristic size of grains for experiments in µm and for simulations in number of cells ETL exp (t), ETL sim (i) equivalent thickness loss as a function of time t for experiments and iterations i for simulations CR exp (t), CR sim (i) corrosion velocity as a function of time t for experiments and iterations i for N X , N Y , N Z dimensions of the system in the CA simulations in number of cells grn ) stationary corrosion velocity of grains as a function of grain corrosion probability λ stat (P grn ) stationary chemical roughness as a function of grain corrosion probability N grn number of grain cells (GRN) in the system N mat number of grain (GRN) and intergranular (IGN) cells in the system N grn-sol number of grain cells in contact with at least a SOL cell N seed number of seeds used to obtain the Voronoï structure

	t, ∆t	time, time duration in experiments in hours h
	i, ∆i	iteration, number of iterations in the simulation
	m(t)	experimental cumulated mass loss as a function of t in g
	ρ	SS density in g.cm -3

simulations v ign (t), v ign (i) grain boundary corrosion rate (corresponding to V gb in article [7]) as a function of time t for experiments and iterations i for simulations v grn (t), v grn (i) grain corrosion rate (V s in article [7]) as a function of time t for experiments and iterations i for simulations α, β groove angles see Figure 4 in degrees • H * groove depth see Figure 4

Table 6 :

 6 List of abbreviations PUREX : Plutonium and Uranium Rening by Extraction

	IGC : intergranular corrosion
	WL: weight-loss
	SROM : stochastic reduced order model
	SS : stainless steel
	AISI : American Iron and Steel Institute
	GB : grain boundaries

HCP : hexagonal close-packed grid CA: cellular automata Cellular automata type of cells: SOL solution, GRN grain, IGN intergranular cells CUDA: Compute Unied Device Architecture

Acknowledgements

The authors would like to thank Damien Féron for his useful comments on this paper, and the CEA and the National Graduate School of Chemistry in Paris for supporting this work.

Data availability

The raw/processed data required to reproduce these ndings cannot be shared at this time as the data also form part of an ongoing study.